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A  major  goal  of  this  program  is  the  selection  and  formulation  of 
suitable  simplified  physics  models.  This  has  been  done  through 
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20.  ABSTRACT  (cont) 

review  of  recent  literature,  discussion  with  mesoscale  Modeling 
experts,  and  a  few  analytic  studies  during  the  first  phase  of  the 
study.  On  the  basis  of  these  studies,  two  different  classes  of 
models  have  been  Identified  for  development  and  implementation. 
A  mass-consistent  wind-field  (or  variational  analysis)  model  has 
been  formulated  that  employs  a  terrain  conformal  coordinate  sya- 
tem.  The  technique  relies  heavily  on  observational  data  and  is 
most  suited  to  a  data  rich  region.  A  computer  program,  VARMET, 
has  been  written,  tests  of  the  method  have  been  performed,  and 
several  comparisons  with  field  data  were  carried  out.  Very  sat¬ 
isfactory  model  performance  was  found  in  these  comparisons.  The 
second  class  of  model  developed  adheres  more  closely  to  first 
principles  and  thereby  is  applicable  to  data-poor  regions.  Based 
on  the  linearized  steady-state  Navler-Stokes  equation,  model 
LINMET  is  written,  and  the  model  performance  is  tested  against 
field  data.  In  addition,  the  primitive  equation  hydrostatic 
code  S1GMET,  which  was  developed  at  SAI  for  application  to  a 
number  of  meteorological  programs,  has  been  adapted  to  the  pres¬ 
ent  program.  Enriched  with  sophisticated  model  physics,  SIGMET 
simulations  form  the  baseline  case  for  two  other  simplified  phys¬ 
ics  models. 

With  cooperation  of  personnel  at  the  US  Army  Atmospheric  Sciences 
Laboratory  on  White  Sands  Missile  Range  (WSMR),  New  Mexico,  data 
were  obtained  from  a  comprehensive  field  observational  program, 
PASS,  which  was  carried  out  at  WSMR  in  late  1974.  These  data, 
together  with  supporting  terrain  height  and  synoptic  weather 
data,  have  been  analyzed  and  Interpolated  for  comparison  with 
corresponding  calculations.  A  second  source  of  comparison  data 
is  the  output  of  calculations  that  take  into  account  more  phys¬ 
ical  phenomena  and  which,  hopefully,  are  in  substantial  agreement 
with  field  data.  Limited  meteorological  and  terrain  data  from 
Fulda,  West  Germany,  were  also  obtained  through  WSMR  and  are  used 
in  a  few  simulations. 
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1 .  INTRODUCTION 


Meteorological  forecasting  of  mesoscale  modifications  of 
synoptic  weather  regimes  is  still  far  from  operational  status. 
However,  within  this  domain  (by  "mesoscale"  we  refer  to  distance 
scales  of  20.-200.  km)  takes  place  most  of  human  activities.  The 
weather  that  people  experience  is  affected  by  their  specific 
location:  within  reach  of  the  land-sea  breeze,  sheltered  or 
exposed  by  terrain,  subject  to  katabatic  winds,  etc.  Consequently, 
the  ability  to  perform  accurate  simulations  on  this  distance  scale 
would  be  highly  beneficial  to  many  activities  that  are  affected 
by  weather.  These  include  the  dispersion  of  pollutants,  transpor¬ 
tation,  siting  of  facilities,  and  the  modification  of  weather 
through  the  inadvertent  or  purposeful  intervention  of  man.  Military 
operations  are  affected  by  most  of  these  weather-dependent  processes 
also.  In  addition  are  the  effects  of  winds  on  the  dispersion  of 
fallout  and  CB  agents,  and  on  the  ballistics  of  projectiles. 

Finally,  military  installations  in  day-to-day  operations  are  deal¬ 
ing  with  the  many  pollution  and  energy  conservation  problems  that 
are  specific  to  each  site.  Mesoscale  weather  affects  these  through 
pollutant  dispersion,  climatological  factors  in  heating  or  cooling 
of  buildings,  and  evaluation  of  renewable  energy  sources  such  as 
wind  and  sun. 

As  the  applications  discussed  above  indicate,  many  (perhaps 
most)  of  the  uses  of  the  mesoscale  simulation  involve  regions  of 
complex  terrain.  While  there  is  a  substantial  history  of  mesoscale 
modeling  research,  this  topic  has  been  neglected  until  recently. 
However,  recently  there  have  been  several  investigations  in  these 
fields.  Aided  by  the  growth  in  computer  capability,  mesoscale 
models  are  now  approaching  the  prototype  application  stage. 
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To  make  further  progress  toward  usable  simulation  tools, 
it  is  now  necessary  to  consider  several  problems  affecting  the 
accuracy  and  economy  of  mesoscale  computer  models.  This  is  being 
done  under  Contract  DAEA18-73-A-0127/0013  entitled  "Improved 
Simulations  of  Mesoscale  Meteorology#"  the  first  phase  of  which 
was  performed  during  the  period  March-October  1977,  and  the  second 
phase  was  performed  during  March-December  1978. 

The  research  program  has  the  objective  of  advancing  the 
technology  of  performing  accurate  and  verified  mesoscale  calcu¬ 
lations  at  a  reasonable  cost  in  computer  usage.  To  do  so,  we  have 
reviewed  the  recent  literature  extensively.  As  an  advanced  start¬ 
ing  point,  the  SAI  SIGMET  primitive  equation  code  has  been  used 
as  the  benchmark  tool  for  mesoscale  simulation.  SIGMET  is  an 
advanced  and  sophisticated  meteorological  mesoscale  code  based  on 
three-dimensional  hydrostatic  equations  in  primitive  form.  The 
mesoscale  experiments  performed  by  the  Atmospheric  Sciences  Labora¬ 
tory  of  WSMR  constitute  a  valuable  data  base  for  the  verification 
of  model  calculations.  These  data  have  been  obtained  and  put  to  use. 
A  search  was  made  for  simpler  computer  codes  applicable  to  flow 
over  complex  terrain,  which  have  greatly  reduced  computing  require¬ 
ments.  Two  such  codes  have  been  selected.  One,  requiring  consid¬ 
erable  local  data  for  initialization,  has  been  implemented  and 
tested.  The  second  formulation,  using  less  data  and  having  more 
physical  content,  is  developed  during  the  second  phase. 

1.1  ROLE  OF  SIGMET  IN  MESOSCALE  MODELING 

The  SIGMET  code  will  be  used  in  two  major  capacities. 

First,  the  code  will  be  applied  to  mesoscale  flow  over  complex 
terrain.  Calculations  will  be  compared  with  field  data  to  indi¬ 
cate  how  accurately  these  simulations  can  be  made.  Second,  the 
SIGMET  code  will  be  compared  with  simplified  physics  code  in  an 
effort  to  reduce  the  expense  of  mesoscale  simulations.  The  com¬ 
parisons  will  permit  trade-offs  to  be  made  between  computing 
speed  and  simulation  accuracy.  The  comparisons  will  also  provide 
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data  for  the  "tuning"  of  the  simpler  codes  through  the  adjustment 
of  parameters. 

1.2  SIMPLE  PRESCRIPTIONS  FOR  OROGRAPHICALLY  INDUCED  FLOWS 

The  presence  of  terrain  in  a  region  of  meteorological 
interest  may  have  large  and  complicated  effects  on  the  local 
wind  field.  The  regional  meteorological  forecaster  is  presented 
with  the  problem  of  inferring  the  effects  of  the  terrain  from  a 
limited  number  of  observations  in  the  vicinity,  from  synoptic 
data,  and  from  previous  experience  in  the  region.  In  the  case 
of  military  operations,  the  amount  of  information  and  experience 
may  be  markedly  less . 

For  almost  every  region  of  interest,  however,  in  contrast 
to  the  sparseness  of  meteorological  data,  there  are  data  on  terrain 
height  in  great  detail.  These  available  data  are  the  basic  ingre¬ 
dient  in  improving  the  objective  analysis  of  the  meteorological 
observations  over  complex  terrain.  The  location  of  mountains  near 
an  observation  site  should  be  taken  into  account  when  those  data 
are  used,  and  perhaps  more  importantly,  mountains  between  observa¬ 
tion  points  should  be  permitted  to  influence  the  estimated  wind 
field  in  their  vicinity. 

As  we  shall  consider  in  the  following  discussion,  several 
reduced-physics  models  of  wind  perturbations  by  mountains  have 
been  developed.  If  these  simple  models  are  able  to  provide  suffi¬ 
ciently  accurate  representations  of  the  perturbation  to  the  wii.d 
field  attributable  to  the  complex  terrain,  a  correction  or  con¬ 
straint  can  be  imposed  on  the  analysis  of  the  observational  data. 

In  the  following  report,  we  describe  in  some  detail  the 
work  that  has  been  performed  under  the  contract.  Data  obtained 
from  the  Atmospheric  Sciences  Laboratory  of  WSMR  have  been  pro¬ 
cessed  for  use  in  comparisons  and  for  -calibration;  this  work  is 
reported  in  Section  2.  The  primitive  equation  code,  SIGMET,  was 
adapted  to  the  investigation  and  calculations  were  performed  with 
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it  as  described  in  Section  3.  In  Section  4,  a  new  3-D  linear  non¬ 
hydrostatic  model  is  described,  and  some  simulation  results  using 
the  new  code,  LINMET,  are  presented.  In  Section  5,  we  report  on 
a  simplified  physics  model  of  flow  over  complex  terrain.  The  char¬ 
acteristics  of  the  simplified  physics  computer  code,  VARMET,  derived 
from  this  study,  are  reported  in  this  section  together  with  exten¬ 
sive  simulations  using  this  model  over  the  White  Sands  region. 
Section  6,  summarizes  the  findings  of  the  study  along  with  some 
future  work.  Documentations  on  the  codes  SIGMET,  VARMET  and  LINMET, 
providing  details  on  code  organization,  subroutine  description, 
input  requirement,  and  sample  setups  are  given  in  Appendix  A,  B,  and 
C,  respectively. 


2.  DATA  ACQUISITION  AND  REDUCTION 


A  major  task  of  the  program  is  to  validate  the  computer 
models  using  available  field  data  and  to  simulate  the  detailed 
wind  field  using  limited  field  data.  Both  the  validation  and  the 
simulation  work  require  intensive  field  data  having  assured  quality. 
The  Atmospheric  Science  Laboratory  at  White  Sands  Missile  Range 
has  provided  us  such  data,  which  are  the  result  of  the  PASS  experi¬ 
ment  conducted  during  November-December  1974.  During  this  intensive 
measurement  program,  rawinsonde  data  were  gathered  at  nine  different 
ground  stations  in  the  region  of  White  Sands  Missile  Range  over 
a  period  of  20  days.  The  reduced  data  were  made  available  in  a 
readily  usable  form  through  magnetic  tape.  In  addition,  terra' 
data  were  provided  for  the  region  on  a  41  x  51  grid  of  points  -h 
5  km  grid  interval. 

2.1  PASS  DATA  OVER  WHITE  SANDS  REGION 

Figures  2.1(a)  to  (d)  show  the  four  prospective  views 
of  the  WSMR  region  using  the  terrain  data.  Figure  2.2  shows  the 
contour  plot  of  the  same.  The  graphs  display  the  Tularosa  Valley 
in  the  center  with  the  San  Andres  Mountains  in  the  west  and  the 
Sacramento  Mountains  in  the  east,  both  running  approximately 
north-south.  These  two  mountain  ranges  have  peaks  up  to  2600  m 
in  height  with  the  valley  floor  averaging  about  1200  m  above  sea 
level.  Figure  2.3  shows  the  ground  locations  of  the  nine  measuring 
stations,  and  Table  2.1  shows  the  UIM  coordinates  and  station  eleva¬ 
tion  (MSL)  for  all  the  stations.  The  meteorological  data  include 
the  pressure,  temperature,  moisture,  magnitude  of  the  horizontal 
velocity  components,  and  the  location  of  the  drifting  balloon  with 
respect  to  the  release  point.  Table  2.2  shows  these  values  for 
a  typical  station  TSX  taken  at/after  0515  MST  on  19  November  1974. 

These  data  have  been  suitably  interpolated  both  vertically 
and  horizontally  to  the  grid  points  used  by  the  finite  difference 


Figure  2.1 


Prospective  Views  of  the  Terrain  in  the 
White  Sands  Missile  Rar.^e  Region.  The 
South-West  Corner  has  the  UTM  Coordinates 
(280  km,  3510  km). 
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Contour  Plot  of  White  Sands  Terrain 
on  a  41x51  Grid  of  Points  and  5  km 
Grid  Interval. 


Figure  2.3  Data  Gathering  Stations  During  PASS  Experiment  in  White  Sands  Region 
(For  names  and  exact  location,  see  Table  2.1.) 


Table  2.1 


Locations  of  Rawinsonde  Stations  Participating 
in  PASS  Experiment . 


Station  Name 

UTM  x-coord 
(km) 

UTM  y-coord 
(km) 

Elevation  MSL 
(m) 

TSX 

375.57 

3586.84 

1232.0 

0R0 

392.11 

3586.57 

1274.0 

MCG 

387.76 

3571.64 

1247.0 

WAR 

366.96 

3572.82 

1216.0 

LSC 

319.80 

3572.70 

1347.0 

SMR 

366 . 36 

3594.73 

1218.0 

RAM 

390 . 66 

3597.13 

1230.0 

APA 

369 . 57 

3610.63 

1204.0 

HMS 

397.69 

3635.88 

1247.0 
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model  for  initialization.  We  have  made  use  of  simple  linear 
interpolation  in  the  vertical  and  !/r2  weighting  in  the  horizontal 
(r  is  the  radial  distance  from  the  point  in  question  to  the  data 
location).  Interpolation  is  done  in  the  transformed  vertical 
coordinate  a  ,  which  is  defined  by 


a 


(z) 


zT-z 


s 


where  z ^  is  the  constant  elevation  of  the  top  surface  of  the 

model,  and  z  is  the  elevation  of  the  terrain.  While  such  a 
s 

transformation  is  appropriate  for  the  simplified  time-independent 
models,  an  alternate  defination  of  o  is  more  suitable  for  the 
more  sophisticated  dynamic  model : 


_  P-PT 

°(P)  =PS-PT 


Here,  p 


T  is  the  constant  pressure  at  the  top  of  the  model,  and 


pg  is  the  surface  pressure.  In  both  definitions  0<o<l  . 

Table  2.3  shows  the  same  measurements  as  Table  2.2,  but  has  two 


additional  entries  a, 


and 


with  z_  chosen  to  be 


(z)  (p)  ’  “T 

6152  m  and  p^  =  500  mbar.  Since  the  balloon  drifts  laterally 

while  rising,  the  appropriate  z  or  p  ,  which  are  the  terrain 

S  S 

height  and  the  surface  pressure  vertically  below  the  balloon 
position,  have  to  be  interpolated.  Since  the  lateral  drift  with 
respect  to  the  release  point  is  known,  this  interpolation  is 
readily  carried  out. 


The  next  step  is  to  construct  tables  similar  to  Table  2.3 
but  containing  values  interpolated  to  the  same  a  -  levels  as 
used  by  the  models.  We  use  the  following  log-linear  formula  for 
generating  a  in  the  models,  which  results  in  fine  vertical  reso¬ 
lution  near  the  surface  and  coarse  resolution  away  from  the  surface. 


12 


£  =  (l-o)  +  a2  £n( 1-o+aj ) 

Here,  values  for  o  are  calculated  corresponding  to  equal 
increments  in  £  .  The  constant  ai  controls  the  first 
increment  in  a  near  the  surface,  and  a?  controls  the  log 
domain.  If  there  are  N  grid  intervals  (N+l  grids)  in  the  ver¬ 
tical,  then  the  increment  in  E,  will  be 

aO)-c(l)  _  l+a2inl(lt_aiiiaii 
A£; - -  N 

Tables  2.4  and  2.5  show  the  same  data  as  Table  2.1  but  interpo¬ 
lated  to  the  a  levels  used  by  the  models.  Table  2.4  is  to  be 
used  by  the  objective  model  and  Table  2 . 5  by  the  predictive  model. 
In  these  tables,  15  vertical' levels  are  generated, and  ax  ,  a2  are 
set  at  0.01  and  0.20,  respectively. 

The  last  step  in  the  procedure  is  to  make  a  horizontal 
interpolation  in  a  -  surfaces  to  initialize  the  values  at  the 
grid  points.  This  is  done  by  utilizing  the  !/r2  weighting 
formula,  which  assigns  maximum  relative  weight  to  the  data  nearest 
to  the  grid  point  and  diminishing  weights  to  the  ones  farther 
away.  In  order  to  initialize  the  meteorological  variable  v  at 
grid  point  (x^^,  y^.^,  °k)  using  data  from  m  stations,  we 
therefore  use 


v 


ijk 


vnk/r 


jjk.nk 


Xw  ^/rijk,nk 
n=l 


where  rijk  nk  is  *he  radial  distance  between  the  grid  point  and 
data  point  in  the  same  a  -  level.  This  is  readily  done  since 
the  UTM  coordinates  of  the  ground  stations,  coordinates  of  the 
data  points  with  respect  to  the  ground  stations,  and  the  coordinates 
of  the  computing  mesh  with  respect  to  the  stations  are  all  known. 
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It  should  be  mentioned  that  there  is  a  time-lag  between 
the  measurements  at  different  levels,  which  is  inherent  in  rawin- 
sonde  type  of  observations.  This  lag  is  about  10  to  15  min 
between  the  ground  measurements  and  measurements  at  500  mbar. 
Moreover,  during  the  experiment  the  data  were  gathered  at  the  nine 
stations  in  two  passes:  four  stations  participating  in  one  pass 
and  the  remaining  five  participating  in  the  second  pass,  a  half- 
hour  later.  Therefore,  there  is  an  additional  30  min  time-lag 
when  one  considers  the  whole  data  set.  In  the  model  initialization, 
we  have  ignored  these  time-lags  and  considered  the  entire  data  as 
representative  of  an  average  state  during  a  1-hour  period. 

Figure  2.4  through  2.11  show  the  velocity  and  temperature 
profiles  chosen  from  the  PASS  experiment  data  for  the  simulations 
over  the  White  Sands  region.  These  data  correspond  to  stations 
TSX ,  ORO,  WAR,  AND  SMR  taken  at  0515,  and  stations  LSC,  RAM, 

APA,  and  HMS  taken  at  0545  on  19  November  1974.  This  day  was 
chosen  from  the  20  days  for  which  measurements  are  available, 
because  the  NWS  weather  record  at  El  Paso,  Texas,  showed  that  the 
day  was  characterized  by  high  wind,  very  little  cloud  cover,  and 
near  neutral  lapse  rate. 

2.2  WEST  GERMAN  DATA 

During  Phase  II,  additional  data  were  made  available 
by  ASL.  These  included  terrain  data  over  a  region  of  West 
Germany  in  Hessen  State  (bordering  East  Germany),  and  upper 
air  data  from  Meiningen,  East  Germany  (east  of  the  above  region). 

The  terrain  data  are  on  a  55  x  40  grid  with  1  km  grid  interval 
with  the  southwest  corner  at  (516.2  km,  5584.3  km)  UTM  coordinate 
with  respect  to  9°  East  Meridian.  Figures  2.12(a)  to  (d)  show 
the  four  perspective  views  of  this  region.  Figure  2.13  shows 
the  contour  plot  of  the  terrain  data.  The  graphs  display  the 
mountainous  and  hilly  terrain  around  Fulda  in  the  state  of  Hessen, 
which  is  typical  of  the  Central  Upland  District  of  West  Germany. 
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Figure  2.4  WSMR  PASS  Experiment  Upper  Air  Measurements,  11/19/74,  Station  TSX. 
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Figure  2.5  WSMR  PASS  Experiment  Upper  Air  Measurements,  11/19/74,  Station  0R0. 
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Figure  2.7  WSMR  PASS  Experiment  Upper  Air  Measurements,  11/19/74,  Station  SMR 
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Figure  2.9  WSMR  PASS  Experiment  Upper  Air  Measurements,  11/19/74,  Station  RAM 
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Figure  2.10  WSMR  PASS  Experiment  Upper  Air  Measurements,  11/19/74,  Station  APA . 
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Figure  2.11  WSMR  PASS  Experiment  Upper  Air  Measurements,  11/19/74,  Station  HMS. 


Figure  2.12  Perspective  Views  of  the  Terrain  Containing 
Fulda,  West  Germany.  The  Southwest  Corner 
has  the  UTM  Coordinates  (516.2  km,  5584.3  km) 
with  Respect  to  9°  East  Meridian. 
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These  mountains  are  small  with  mild  slopes.  The  highest  mountain 
in  this  region  is  Wasserkuppe  Mountain  on  the  southeast  boundary, 
with  its  peak  elevation  at  950  m  above  sea  level.  The  Vogels- 
terg  mountain  range  is  in  the  southwest  area  of  this  region, 
with  its  highest  peak  at  774  m,  which  is  outside  the  region  con¬ 
sidered.  The  city  of  Fulda  is  in  a  valley  located  near  the 
southeast  part  of  this  region,  with  the  valley  floor  elevation 
ranging  from  250  m  to  400  m.  The  Fulda  River  streams  past  the 
city  to  the  north,  separating  the  mountainous  regions. 

For  simulation  over  this  region,  as  also  with  the  White 
Sands  region,  a  limited  area  containing  32  x  32  grid  points  is 
chosen.  For  the  West  Germany  region,  an  area  with  its  southwest 
corner  at  (526.2  km,  5584.3  km),  UTM  coordinate  is  taken.  Since 
simulations  over  this  region  were  carried  out  using  the  linear 
model  LINMET  which  presumes  periodic  boundary  conditions  in  both 
lateral  directions,  a  boundary  smoother  is  applied  over  a  5-grid- 
wide  strip  all  around  the  limited  area  in  order  to  avoid  abrupt 
slopes  at  the  boundaries.  This  smoother  sets  the  boundary  eleva¬ 
tion  at  a  constant  value  of  400  m  and  then  linearly  interpo¬ 
lates  the  elevation  between  the  first  and  the  fifth  grid.  The 
elevations  at  the  fifth  grid  and  beyond  are  left  at  the  original 
values.  This  smoother  effectively  removes  the  boundary  problem 
and  seems  to  work  satisfactorily.  Figures  2.14(a)  and  (b)  show 
the  perspective  and  contour  view  of  this  limited  area.  Figures 
2.15(a)  and  (b)  show  the  same  for  the  limited  area  in  the  White 
Sands  region  over  which  some  LINMET  simulations  are  carried  out. 
For  this  region,  the  boundary  elevation  is  maintained  at  1200  m 
above  sea  level,  and  the  southwest  corner  is  located  at  (340  km, 
3570  km)  UTM  coordinates.  Figures  2.16  and  2.17  show  two  upper 
air  measurements  taken  at  Meiningen,  East  Germany,  which  are  used 
for  a  sample  simulation  over  the  Fulda  region.  Meiningen  is 
located  east  of  the  Fulda  region. 


29 


Figure  2. 14 


Perspective  (a)  and  Contour  (b)  Views  of  the 
Region  Around  Fulda,  West  Germany,  on  a  32  x  32 
Grid  with  1  km  Grid  Interval.  The  Southwest 
Corner  has  UTM  Coordinates  of  (526.2  km,  5584.3 
km)  and  the  Boundaries  are  Smoothed. 


Figure  2.15  Perspective  (a)  and  Contour  (b)  Views  of  the 

White  Sands  Region  on  a  32  x  32  Grid  with  5  km 
Grid  Interval .  The  Southwest  Corner  has  UTM 
Coordinates  of  (340  km,  2570  km)  and  the 
Boundaries  are  Smoothed. 
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Measurement  at  Meiningen,  East  Germany. 


Figure  2.17  Upper  Air  Measurement  at  Meiningen,  East  Germany 
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o.  PRIMITIVE  EQUATION  MODEL  SIGMET 


This  class  of  wind-field  calculation  is  defined  as  solv¬ 
ing  the  primitive  equations.  Similar  systems  of  equations  are 
employed  for  calculations  of  global  circulation,  synoptic  weather 
forecasting,  mesoscale  forecasting,  and  for  research  in  severe 
storms  and  cumulus  convection  modeling.  in  this  formulation,  it  is 
quite  straightforward  to  include  the  effects  of  boundary  layer 
turbulence  and  thermal  forcing  (as  is  done  in  the  SIGMET  computer 
code  of  SAI).  These  models  are  applied  by  performing  a  3-D  time 
dependent  calculation  starting  from  initial  conditions  and  employ¬ 
ing  boundary  conditions  determined  by  the  synoptic  or  larger-scaie 
weather  regime.  In  the  case  of  diurnally  cyclic  flow,  the  calcu¬ 
lation  should  continue  for  a  day  or  more  to  determine  the  distri¬ 
bution  of  winds  at  each  site.  Two  different  treatments  of  the 
vertical  velocity  can  be  distinguished,  based  on  the  magnitude  of 
the  expected  vertical  acceleration  in  the  flow  field.  If  the 
vertical  acceleration  is  small,  as  is  the  case  when  the  size  of 
zones  in  the  horizontal  direction  is  somewhat  greater  than  1  km 
and  when  the  terrain  is  relatively  smooth,  the  hydrostatic  approx¬ 
imation  can  be  employed.  This  simplifies  the  calculation  and 
makes  the  transformation  to  a  terrain-conformal  coordinate  system 
(as  used  in  SIGMET)  easier.  If  the  horizontal  spacing  is  small, 
the  vertical  acceleration  equation  can  be  retained.  Frequently 
the  resulting  equations  are  solved  in  the  Boussinesq  approximation, 
in  vhich  the  effects  of  buoyancy  are  present  only  in  the  equation 
for  the  vertical  component  of  momentum. 

The  primitive  equation  models  are  the  only  class  in  which 
all  of  the  physical  processes  affecting  near  surface  winds  are 
included.  Other  considerations  aside,  the  primitive  equation 
models  are  the  obvious  choice  for  mesoscale  applications.  Factors 
weighing  against  their  use  are:  complexity,  incomplete  verifica¬ 
tions  of  the  codes,  and  expense.  Presently,  the  calculation  of 


one  case  of  a  3-D  mesoscale  wind  field  costs  a  few  hundred  dollars 
for  a  short  calculation,  but  this  cost  will  be  reduced  by  code 
optimization,  faster  computers,  and  the  technique  of  partial  impli¬ 
es  tation.  The  combination  of  these  could  reduce  the  cost  by  a  fac¬ 
tor  of  10  or  more. 


3.1  SIGMET  FORMULATION 

The  SIGMET  family  of  mesoscale  (10  ->100  km)  meteorology 
simulation  codes  is  based  on  the  work  of  Anthes  (1972)  and  uses 
a  finite  difference  technique  to  solve  the  so-called  primitive 
equations  t'  describe  transient  flow  in  the  atmosphere.  One-  two- 
and  three-dimensional  versions  of  the  code  have  been  written,  each 
of  which  share  the  same  level  of  physical  approximation  as  well  as 
similar  numerical  techniques.  These  common  features  among  the 
codes  have  permitted  the  development  and  testing  of  much  of  the 
physical  and  numerical  modeling  using  the  more  economical  1-D  or 
2-D  versions  with  subsequent  implementation  in  the  3-D  code.  Also, 
although  the  3-D  version  is  required  ultimately  for  application  to 
mesoscale  simulations  of  regional  meteorology,  it  is  believed  that 
much  can  be  learned  with  the  simpler  and  more  economical  codes. 

The  SIGMET  codes  solve  the  conservation  equations  for  the 
relevant  meteorological  variables  (wind  components,  temperature, 
and  moisture)  and  describe  both  the  planetary  boundary  layer  and 
upper  level  flow  as  they  are  affected  by  complex  terrain.  The 
equations,  formulated  in  an  Eulerian  framework,  account  for  advec- 
tion,  Coriolis  effects,  turbulent  heat,  momentum  and  moisture  trans¬ 
port,  and  viscosity.  In  many  of  the  applications  of  mesoscale 
simulation,  the  atmospheric  boundary  layer  winds  are  of  particular 
interest.  Over  regions  of  complex  terrain,  these  low  level  winds 
suffer  from  numerical  errors  when  finite  differences  are  formed 
in  Cartesian  Coordinates.  In  order  to  avoid  this  problem,  a 
coordinate  transformation  has  been  carVied  out  to  a  new  coordinate 
system  (the  sigma,  a,  coordinate  system)  in  which  the  lowest 
coordinate  surface  is  conformal  to  the  terrain  surface.  In  this 
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representation,  after  care  is  exercised  in  calculating  the  trans¬ 
formed  pressure  gradients,  the  needed  resolution  and  accuracy  near 
the  surface  are  achieved  by  concentrating  the  finite  difference 
grid  in  the  near  surface  region.  Expanding  the  grid  into  the 
upper  levels  with  a  non-constant  grid  distribution  satisfies  the 
less  stringent  accuracy  requirements  there  while  minimizing  ver¬ 
tical  grid  points  to  maintain  an  economical  calculation. 

Diurnally  varying  winds  such  as  the  land-sea  breeze, 
slope  winds,  etc.,  are  important  in  the  regional  meteorology. 

These  winds  are  caused  by  thermal  forcing  of  the  mesoscale  flow, 
and  the  simulation  requires  that  solar  radiation  and  its  related 
phenomena  be  included  in  the  model.  To  do  so,  subroutines  of 
SIGMET  have  been  developed  in  which  solar  and  terrestrial  radiation 
are  calculated.  These  affect  the  winds  by  inducing  thermal 
pressure  gradients  across  regions  of  differing  albedo  and  along 
sloping  terrain.  These  radiation  terms  are  also  affected  by  the 
amount  of  atmospheric  water  vapor  and  by  the  conduction  and 
storage  of  heat  in  the  soil.  Consequently,  soil  temperature  and 
moisture,  and  atmospheric  water  vapor  mixing  ratio,  which  satisfy 
prognostic  equations, are  included  in  SIGMET.  The  development  and 
testing  of  the  surface  heat  and  mass  transfer  balances  have  resulted 
in  significant  improvements  in  the  physical  accuracy  and  numerical 
economy  of  their  treatment.  As  a  result  of  recent  work,  all  ver¬ 
sions  of  the  SIGMET  code  are  currently  at  the  same  level  of  devel¬ 
opment  in  terms  of  the  physical  modeling  items  mentioned. 

A  very  important  physical  effect,  which  must  necessarily 
be  taken  into  account  in  mesoscale  calculat ions, is  the  transport 
of  momentum,  energy,  and  water  vapor  by  turbulence  in  the  atmo¬ 
spheric  boundary  layer.  This  is  a  field  of  active  current  research 
in  fluid  dynamics.  Several  different  approaches  to  the  description 
of  turbulent  transport  have  been  studied,  and  the  best  compromise 
for  meteorological  flows  seems  to  be  the  "level  2"  and  "level  2£" 
models  of  Yamada  and  Mellor  (1975).  The  algebraic  equations 
designated  as  a  "level  2"  model  by  the  authors  correspond  to  a  mixing 
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length,  equilibrium  turbulence  theory  in  which  full  account  of 
density  stratification  effects  is  retained.  Comparison  of  the 
level  2  predictions  of  the  atmospheric  boundary  layer  with  those 
of  more  complicated  approximations  indicates  that  they  capture  many 
of  the  features  of  boundary  layer  turbulence  (Mellor  and  Yamada, 
1974).  This  purely  local  theory,  however,  does  not  account  for 
turbulence  history  effects,  i.e.,  nonequilibrium  effects  caused 
by  flow  over  terrain,  changes  in  surface  roughness,  etc.  Such 
effects  are  contained  in  the  level  2§  model  which  as  used  in  the 
present  work  includes  a  turbulent  energy  transport  equation. 

Because  of  the  importance  of  such  effects  in  complex  terrain  to 
the  structure  of  the  boundary  layer  winds,  the  level  2$  model  has 
recently  been  studied  and  incorporated  in  all  versions  of  the 
SIGMET  code.  Barring  unforeseen  difficulties  with  the  model,  this 
completes  the  turbulence  modeling  work, since  it  is  believed  that 
the  present  model  is  as  detailed  as  required  in  the  present  program. 

Strong  turbulent  diffusion  across  the  small  zones  required 
to  resolve  the  planetary  boundary  layer  near  the  surface  may,  under 
some  circumstances,  limit  the  time  interval  to  a  few  seconds  if  an 
explicit  formulation  is  used.  In  order  to  remove  this  limitation, 
a  new  formulation  of  the  vertical  diffusion  which  is  partially 
implicit  has  been  implemented.  This  method  eliminates  the  time 
interval  restriction  associated  with  diffusion.  As  indicated  by 
a  linear  stability  analysis,  it  also  markedly  relaxes  the  stability 
inequality  associated  with  vertical  advection.  The  partially 
implicit  code  contains  simultaneous  linear  equations  in  each  verti¬ 
cal  column  of  zones.  These  equations  are  equivalent  to  a  coeffi¬ 
cient  matrix  of  the  unknown  quantities  which  is  tridiagonal  and 
is  solved  by  a  forward/backward  substitution  algorithm.  This  ver¬ 
tical  diffusion  subroutine  has  been  incorporated  into  all  versions 
of  SIGMET,  and  the  resulting  code  has  been  debugged  and  tested. 

In  summary,  the  SIGMET  modeling ’contains  all  of  the  physics 
required  to  simulate  the  time  dependent  meteorology  of  a  region 
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based  upon  minimal  weather  station  data.  However,  the  relative 
expense  of  the  calculations,  together  with  desire  to  increase  the 
amount  of  wind-field  modeling  possible  within  a  limited  budget, 
has  prompted  the  continuing  efforts  to  increase  the  numerical 
efficiency  of  the  codes. 

We  now  outline  the  equations  that  are  solved  in  the 
3-D  SIGMET  code.  These  consist  of  the  equations  for  the  3-D 
distributions  of  the  horizontal  components  of  the  wind  velocity, 
the  potential  temperature, and  the  mixing  ratio  of  the  water 
vapor.  In  addition,  there  is  an  equation  for  the  2-D  distribution 
of  the  pressure  thickness  of  the  atmosphere.  In  terms  of  these 
quantities  the  equivalent  vertical  velocity  a  is  obtained  as  a 
diagnostic  quantity. 

By  virtue  of  the  condition  of  hydrostatic  balance,  it  is 
possible  to  simplify  the  calculation  substantially  by  eliminating 
the  differential  equation  for  the  vertical  velocity.  In  addition, 
it  proves  useful  to  transform  the  vertical  coordinate,  replacing 
the  altitude  with  a  scaled  pressure  coordinate.  This  pressure 
coordinate,  called  the  sigma  coordinate,  is  defined  as  follows: 


P  -  P„ 


The  hydrostatic  relation 


(3.1) 


dz  _  _  _n_ 
da  gp 


(3.2) 


gives  the  relationship  between  the  altitude  and  the  a  coordinate 
The  sigma  coordinate  is  in  the  range  0  s  o  5  1  ,  and  takes  the 
value  a=l  at  the  surface.  It  is  applicable  to  a  surface  with 
complex  terrain,  as  illustrated  in  Figure  3.1. 
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(a).  Two-Dimensional  Terrain  in  x,z  Space. 


(b).  Two-Dimensional  Terrain  In  x,z  Space  With 
Contours  of  Constant  e  Superimposed. 


(e).  Two-Dimensional  Terrain  in  x,a  Space.  The 

Location  of  the  Terrain  is  Indicated  by  the  Cross- 
Hatched  Region. 


Figure  3.1  Terrain  Coordinate  Transformation. 


This  transformation  is  used  to  eliminate  the  altitude. 
After  considerable  manipulation,  the  primitive  equations  take 
the  form 
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(3.3) 


The  effects  of  vertical  and  horizontal  turbulent  diffusion 
are  contained  in  the  terms  Fy  and  FH  of  Eq.  (3.3).  These  are 
approximated  by  Fickian  diffusion  expressions. 
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Eqs.  (3.4)  and  (3.5)  contain  diffusion  coefficients  for 
vertical  and  horizontal  transport  by  turbulent  eddies.  The  hori¬ 
zontal  coefficients  have  been  evaluated  from  empirical  data,  and 
the  vertical  diffusion  coefficients  are  evaluated  using  Yamada- 
Mellor  turbulence  closure  theories. 
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The  first  two  equations  of  Eq .  (3.3),  describing  the  rate 
of  change  of  the  horizontal  velocity  components,  contain  terms  for 
the  horizontal  pressure  gradient  and  the  Coriolis  force.  Each 
component  of  the  pressure  gradient  consists  of  two  terms:  the 
first  is  proportional  to  the  gradient  of  geopotential  height,  and 
the  second  contains  the  gradient  of  the  surface  pressure.  These 
two  terms  are  evaluated  by  independent  difference  approximations 
which,  over  complex  terrain,  will  contain  truncation  errors.  We 
have  found  that  these  errors  can  be  appreciable  and  have  devised 
a  scheme  for  reducing  them  to  a  negligible  value.  The  Coriolis 
parameter  f  =  2fi  sine  can  be  taken  to  be  constant  over  the  meso- 
scale  region.  The  third  of  Eq.  (3.3)  is  for  the  rate  of  change  of 
the  energy.  In  Q,  we  take  into  account  at  present  only  the  radia¬ 
tive  flux  divergence,  neglecting  the  liquid  water  effects. 

Difference  equations  approximating  the  above  differential 
equations  attempt,  so  far  as  possible,  to  retain  second-order 
accuracy  in  the  spatial  and  time  dependence.  The  spatial  differ¬ 
ence  equations  employ  a  rectangular  mesh  in  which  velocity  vari¬ 
ables  are  offset  from  the  thermodynamic  variables.  The  time 
derivatives  of  the  equations  employ  three  levels  in  the  finite 
difference  scheme.  This  is  basically  a  leapfrog  formulation. 


3.2  SIGMET  SIMULATIONS 

Relatively  short  duration  simulations  of  the  wind  field 
around  White  Sands  Missile  Range  have  been  performed  using  the 
three-dimensional  version  of  SIGMET.  The  calculations  utilized 
the  limited  terrain  shown  in  Figures  3.2  and  3.3.  This  region 
encompasses  all  the  measurement  locations,  major  parts  of  the 
Tularosa  Valley  and  the  San  Andres  Mountains,  and  parts  of  the 
Sacramento  Mountains.  The  finite  difference  grid  design  is 
defined  in  the  following. 
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This  results  in  a  solution  domain  of  125  km  (east-west)  by  105  km 
(north-south)  around  WSMR,  which  is  defined  by  the  borders  of  the 
contour  plot  mentioned  above.  The  relatively  coarse  5  km  horizontal 
zoning  was  chosen  in  an  effort  to  include  as  much  of  the  wind  field 
influencing  terrain  around  White  Sands  as  possible,  while  retaining 
a  relatively  economical  calculation.  The  simulation  was  initialized 
at  0515  LST  using  the  wind  speed  components,  temperature,  and  water 
vapor  concentration  given  by  the  station  TSX  upper  air  sounding  at 
that  time  on  19  November  1974  (see  Figure  2.4). 


From  this  initial  state,  the  equations  were  integrated 

with  a  time  step  of  3  sec  to  a  total  simulation  time  of  7£  min 

(or  150  cycles).  The  3  sec  time  step  is  considerably  below  the 

maximum  time  step  allowed  by  the  gravity  wave  stability  limitation 

of  the  finite  difference  scheme  in  SIGMET,  i.e.,  At  <  min  (Ax, Ay)/ 

c 

2/gH  ,  where  H  is  the  maximum  depth  of  the  atmospheric  layer  being 
simulated.  In  order  to  test  the  timewise  truncation  error  in  the 
solution  (formally  second-order),  the  calculation  was  repeated  with 
At  =  9  sec ,  which  is  just  under  (80  %  of)  the  critical  time  stepAtc. 
At  a  simulation  time  of  7i  min,  the  resulting  field  variables  from 
each  calculation  compared  to  within  -vl.%  indicating  the  absence  of 
significant  timewise  truncation  error. 


The  results  of  the  At  =  3  sec  calculation  are  presented 
in  Figures  3.4  through  3.7  below.  The  near  surface  wind  vectors 
and  speed  contours  shown  in  Figures  3.4,  for  ^>15. m  AGL,  and 
Figure  3.5,  ^200. m  AGL,  illustrate  the  development  of  the  expected 
terrain  effects.  They  show  the  speed-up  on  the  lee  side  of  the  San 
Andres  Mountains,  near  the  center  of  the  region,  and  the  turning 
and  flow  acceleration  due  to  the  higher  Sacramento  Mountains  to  the 
east  and  northeast.  Also  shown  is  the  tendency  toward  low  winds 
in  the  central  valley.  Finally,  additional  details  of  the  simulation 
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Figure  3.2  Prospective  View  of  Limited  White  Sands  Terrain  (120  km  x  100  km)  on 
5  km  grid  used  for  SIGMET  Simulation.  The  South-West  corner  has  the 
Coordinates  (310  km,  3540  km). 


Figure  3.3  Contour  Plot  of  Limited  White  Sands  Terrain 
Used  for  SIGMET  Simulation.  The  Rectangle 
Shown  by  Dashed  Lines  is  used  for  VARMET 
Simulat ions . 


are  presented  in  the  near  surface  temperature  contours,  Figure 
3.6;  water  vapor,  Figure  3.7;  and  turbulent  energy,  Figure  3.7. 
The  turbulence  energy  contours  show  the  high  turbulence  in  the 
high  shear  regions  of  the  flow;  in  addition,  they  show  the 
relatively  low  values  in  the  low  wind,  stable  valley  floor  areas 
of  the  flow. 


Figure  3.4  Near  Surface  (^15  m  AGL)  Wind  Field 
Vectors  and  Speed  Contours :  SIGMET 
Simulation  of  WSMR,  0515,  11/19/74. 
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Figure  3.7  Near  Surface  Water  Vapor  Mixing  Ratio 

SIGMET  Simulation  of  WSMR ,  0515,  11/19/74 


Figure  3.8  Near  Surface  Mean  Square  Turbulent  Velocity 
Fluctuation  (m2/sec3);  SIGMET  Simulation  of 
WSMR,  0515,  11/19/74. 
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LINEAR  MODEL  LINMET 
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Based  on  steady  linearized,  inviscid  Navier-Stokes 
equations  in  three-dimensions,  a  model  is  developed  to  predict 
wind  field  in  data  poor  regions  with  moderate  success.  While 
this  class  o'  models  will  not  simulate  the  detailed  physics 
of  the  i ‘ow, like  the  sophisticated  primitive  equation  models 
described  in  the  previous  chapter,  the  computational  effort 
will  be  significantly  less.  It  may  thus  be  possible  to  use 
such  models  in  real-time  operations.  One  such  model,  LINMET, 
is  developed  here.  The  theoretical  basis  for  its  formulation 
and  some  sample  calculations,  using  the  code,  are  given  in  the 
following . 


4.1  LINMET  FORMULATION 

In  order  to  construct  a  compact  and  self-consistent  set 
of  finite  difference  equations  for  LINMET,  we  discuss  both  the 
differential  equations  of  motion  and  the  discretizations  of  the 
equations  for  numerical  computation  side  by  side.  Extensive  use 
of  Fast  Fourier  Transform  (FFT)  is  contemplated  in  this  model. 

We  assume  that  the  terrain  height  is  represented  at 

the  Cartesian  Coordinates  xi  =  XQ  +  i^x,  =  yQ  +  jAy  (i  =  1,  N] 
j  =  1,  ...N2)  separated  by  constant  intervals  Ax,  Ay.  The  orien¬ 
tation  of  x  and  y  axes  is  arbitrary.  Assuming  that  the  boundary 
conditions  are  periodic  in  both  x  and  y,  then  the  finite  Fourier 
transform  HkJl  can  be  readily  obtained  from  the  data; 


”k*  =  SIX  6XP  +  (j-l)(*-D/N2]} 

1=1 j=l 


(4.1) 


k  =  1,2  ...  Nj,  i  =  1,2  .*..  N2  and  l  =  /-I 
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The  inverse  transform  is  given  by  a  similar  expression; 


|-2i.-t[(i-l)(k-l)/N1 


+  (j-l)(i-l)/N2]}.(4.2) 


Since  these  finite  representations  do  not  admit  differentiation, 
it  is  necessary  ultimately  to  form  difference  equations  for  the 
atmospheric  variables  that  can  be  represented  in  the  form  indi¬ 
cated  in  Eqs .  (4.1)  and  (4.2).  Consequently,  in  the  following 
discussion,  we  derive  the  equations  satisfied  by  components  of 
the  discrete  Fourier  transform  of  the  fluid  dynamic  variables. 

The  linearized  steady-state  equations  for  the  conservation 
of  mass,  momentum  and  energy  of  the  atmosphere  are: 
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(4.3) 


In  the  above  equations,  differentiation  with  respect  to  the 
horizontal  coordinates  x,  y,  or  the  vertical  coordinate  z,  is 
represented  by  the  corresponding  subscript.  Unperturbed  vari¬ 
ables  (corresponding  to  horizontal,  steady  flow  in  the  absence 
of  terrain)  are  denoted  by  U,V,  and  the  overbarred  quantities 
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that  depend  only  on  z.  The  dependent  variables  are:  u,v,w,  the 
wind  vector  components;  and  p,p,T,  the  air  density,  pressure,  and 
temperature.  The  other  quantities  appearing  in  Eq .  (4.3)  are 
the  adiabatic  lapse  rate  r  =  g/Cp,  the  specific  heat  of  air  at 
constant  pressure  Cp,  and  the  gravitational  constant  g.  The 
unperturbed  atmosphere  is  in  hydrostatic  equilibrium  with  p  - 
-  gp.  The  above  atmospheric  variables  are  related  to  the  terrain 
height  through  the  linearized  boundary  condition 

w(o)  =  U(o)H  +  V(o)H  .  (4.4 

x  y 

In  order  to  reduce  Eq.  (4.3)  to  a  more  compact  form,  we  eliminate 
T  and  p,  obtaining 

wz  +  (S-s)  w  +  ux  +  vy  +  (upx  +  vPy)/( pC2  )  =  0 

U2wxx  +  2UVwxy  +  V2wyy  +  gSw  +  (Upxz  +  Wpyz^*  (4'5 

+  g(Upx  +  VPy)/(pC2)  =  0 

where  S  =  (T  +D/T 
z 

s  =  -Pz/? 

and 

C2  =  (C  p)/(C  p),  the  adiabatic  sound  speed, 
s  p  v 

These  equations,  together  with  the  second  and  third  equations 
of  Eq.  (4.3),  constitute  four  equations  for  the  quantities  u,v,w, 
and  p.  Rather  than  performing  further  reduction  of  the  equation 
in  differential  form,  requiring  additional  differentiation,  we 
proceed  to  form  difference  equations  that  are  subsequently 
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subjected  to  discrete  Fourier  transform.  The  resulting  equations 
may  be  solved  algebraically.  By  applying  Fourier  transform  to 
the  discretized  equations  rather  than  to  the  continuous  equations, 
we  avoid  the  negative  wave  numbers.  Moreover,  it  can  be  shown 
that  this  approach  results  in  more  accurate  representation  of 
higher  wave  number  modes. 

Eqs.  (4.3)  and  (4.4)  contain  first  and  second  derivative 
terms,  which  are  to  be  approximated  by  central  difference  expres¬ 
sions.  We  define  each  of  the  variables  at  the  same  position 
xi ’  y i ‘  Forming  centered  differences,  we  obtain  the  following 
approximations . 


(  ^x) ij 


=  *i+l,j  ~  *1-1,  j 

2Ax 


L  )  =  J _ .t -l,j 

\  xx/ ii  .  , 


AX‘ 


L  )  -  L  )  -  ~  h-l.J+1  '  »i+l.J-l  *  ♦i-l.J-l  .  (1  e) 

'  xy'  ij  '9yx'  ij  4AxA y 


We  now  transform  these  terms  using  Eq .  (4.1)  and  taking  into 
account  that  the  boundary  conditions  are  cyclic.  For  example, 


i  AN2 

(*x)ij~2A^^2  (♦i+l.j  "  ♦i-lj) 


N,  N2 
1=1 j=l 

exp  2^|[(i-l)(k-l)/N1  +  (j-1)  (  £-1  )/N2  ]} 


_  zL  Sin 1 2w(k-l  )/Ni 1  * 

Ax  uk  t 


*  Sk*ke  ‘ 


53 


Similarly,  we  obtain 


(  *y|ij  *  1  Vkl 

(*xx  )ij  "  ^k^ka 
(*yy  )ij  =  Ci*kt 

and  ( **y  lu  =  •  skVk£  • 

Where  Sk  =  -  Sin  [2tt  (k-l)/Ni]/£x 
S4  -  -  Sin  [2rr  ( £-l)/N?]/Ay 

Ck  =  aI'0  {Cos  0(k-l)/Nj]  -l} 

and  =^2  {Cos  [2* (  t-1  )/N2  ]  -l}  . 

Using  Eq.  (4.7),  the  second  and  third  equations  of  Eq.  (4.3) 
and  both  equations  of  Eq.  (4.5)  can  be  transformed  to  the  dis¬ 
cretized  Fourier  space. 

*  ~  Sk  “ 

V  '  *Uzw  +  —  p  =  0 

U  V  -  -tV  W  +  -=r-  P  =  0 
n  z  p 

wz  +  (S-s)w  +  ^Sku  +  SlV+  p)  -  o 

(u2Ck  -  2SkS^UV  +  V2C!L  +  gs)  w  +  >iy  (  Pz  +  c^p)*0 

s 


(4.7) 


(4.8) 
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where  U  *  US.  +  VS. 

n  k  l 

and  the  kt  subscript  has  been  omitted  from  each  of  the  transformed 
variables . 

From  the  first  three  equations  u  and  v  can  be  eliminated  to 

give 


By  elimination  of  p  and  p  between  Eq .  (4.8)  and  (4.9),  we  obtain 

A  Z 

a  single  equation  for  w. 


w 


zz 


-  (•  -  v)  »z  *  («■'*  -  K2) 


w  *  0 


where 


U«j2  '  Vi  /  u s  Vs  +  a  /  U” 


R2  •  (-  “2Ck  +  J\s(m  -  v2°t)  o 


(4.10) 


_  Un2 

sk2+s^ 


a  -  1 


(4.11) 


The  equations  derived  above  are  based  on  a  difference  formulation, 
which,  for  small  wave  numbers,  would  be  expected  to  limit  to  a 
more  familiar  differential  expression.  Consequently,  let  us 


examine  the  behavior  of  the  following  coefficients  when  the 
indices  k  and  l  are  small. 


1 


sin[ 2'i  ( k-1 ) /Ni  ]  _  2Ti(k-l 
Ax  "  N i Ax 


{cos  [2.(k-1)/N1]  -  1}  .  -[2^|^1)]  - 


Similarly , 


(4.12) 


-  k  2 

y 


The  discrete  wave  numbers  k  and  k  are  defined  in  terms  of  the 

x  y 

indices  k  and  l  and  the  lengths  x  =  NjAx  and  y  =  N2Ay.  These 
limiting  values  can  be  introduced  in  Eq.  (4.11). 


k  1)  +  k  V 

V  *  ^-V— V  =  u  : 

k*  +  ky 


■t-=4) 


(4.13) 


2 

K2  *  (kx2»2  *  2kXkyVV  +  VV!)  i1  -  57)  1  V 

l'2  «  (s  +  _  (u.n.),a  +  (s  +  -z)  lHaU 

“n2  V  V  }  Un  \  V  /  U„ 


r 


These  terms  give  the  limiting  form 


,s  +  ~  wz  +  L 


L'2  -  (kx7  *  V)  ('  -$)  »  -  °- 

v  S  /J  (4.14) 


In  order  to  place  Eq.  (4.10)  in  a  form  convenient  for 
numerical  integration,  we  eliminate  the  first  derivative  term. 
Introducing  the  modified  vertical  velocity  term. 


p _  p  ( o) 


0  (  O  )  P 


we  obtain  the  familiar  Scorer  equation  in  canonical  form. 


Wzz  +  [L2  -  K2]  W  =  0 


where 


-  T  '  2  4.  1  fe  +  Z 

-  L  ‘  +  7  Is  +  

4  \  p 


(4.15) 


(4.16) 


In  order  to  find  the  horizontal  components  of  the  perturbed  wind 
velocity,  the  first  two  equations  of  Eq.  (4.8)  and  (4.9)  are  used. 
Substituting  for  pressure,  we  obtain 


V2  ,  L  .  lus)z\ 

“  ~ 1  iwl  ^  w| 

:  t  I  V,  L  .  (usL\  sn  ,  vz|  -) 

Consequently,  these  quantities  can  be  obtained  directly  when 
w  and  its  vertical  derivative  w„  are  known. 

Z 


. 


A 


(4.17) 
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Boundary  conditions  arc  required  at  the  top  and  bottom  of 
the  computational  region.  These  conditions  reflect  the  physical 
constraints  that  are  imposed  at  these  locations.  At  the  surface, 
the  flow  is  assumed  to  follow  the  terrain  without  taking  account 
of  dissipation  in  the  boundary  layer.  In  the  linear  approximation, 
the  boundary  condition  is  applied  at  z  =  0.  Here  the  unperturbed 
horizontal  velocity  is  deflected  by  the  sloping  terrain  inducing 
a  vertical  velocity. 


w(o)  =  U(o )  H  +  V(o)  H,  . 

x  y 


(4.18) 


We  need  the  boundary  condition  for  W(o),  which  is  derived  from 
Eq.  (4.18)  by  replacing  the  derivatives  by  central  difference 
approximations  and  performing  the  discrete  Fourier  transform 


W(o)  =  i.  [U(  o)Sk  +  V(o)S£]HkJl 


=  Hk* 


(4.19) 


At  the  upper  boundary,  the  physical  condition  is  not  understood  as 
well.  It  is  known  that  certain  gravity  waves  may  propagate  into 
the  stratosphere  where  at  high  altitude  ('vlOO  km)  they  are  dis¬ 
sipated.  Under  other  circumstances , the  gravity  waves  may  be  sub¬ 
stantially  or  completely  reflected  by  adverse  wind  gradients. 

These  two  alternatives  are  crudely  incorporated  into  the  formula¬ 
tion  by  evaluating  whether  the  solution  displays  exponential  or 
oscillatory  behavior  at  the  upper  boundary.  In  the  former  case, 
corresponding  to  the  coefficient  of  W  in  Eq .  (4.16) 


(■’  -  '■%; 


>  0 


(4.20) 


The  boundary  conditions  are  given  by 


-  \/k2  - 


L2  W 
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For  the  second  case,  when 


<  0 


the  boundary  conditions  are  chosen  to  correspond  to  upward  propa¬ 
gating  waves  alone.  These  conditions  are  achieved  when 


5 

For  numerical  solution,  the  second-order  ordinary  differ¬ 
ential  Eq .  (4.16)  is  decomposed  into  a  set  of  four  first-order 
equations  (corresponding  to  real  and  imaginary  parts  of  W)  and 
are  solved  using  sixth-order  Adams -Moult on  Predictor-Corrector 
method,  with  starting  procedure  based  on  Roser’s  formula.  The 
Fourier  analysis  and  synthesis  are  performed  by  a  version  of 
the  Cooley-Tukey  algorithm.  The  one  we  use  requires  the  dimen¬ 
sion  of  the  data  in  both  directions  to  be  an  integral  power  of 
two.  While  it  is  possible  to  acquire  a  generalized  algorithm 
for  any  arbitrary  dimension,  the  computational  expense  will  be 
higher . 

In  order  to  solve  the  boundary  value  problem  Eq .  (4.16), 
subject  to  conditions  Eqs.  (4.19)  and  (4.20)  or  (4.21),  it  is 
advantageous  to  reformulate  the  equation  into  an  initial  value 
problem  for  the  purpose  of  numerical  integration.  First  we  notice 
that  a  general  solution  of 

V  +  (l2  -  K2)  i  -  0 
is  given  by 

W  =  aUi  +  bU2 

with  u,  (zT)  -  1,  ulz  (zTj  =  0 


(4.21) 


W  =  i /Z7  -  K2  W  sgn 
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and 


1  . 


(4.22) 


The  quantity  W  is  complex  and  both  real  part  WR  and  imaginary 
part  Wj  obey  the  same  equation  (since  L?  -  K2  is  real)  and  are 
given  by 


R 

I 


aRul 


bRU2 


aT  u 
I  1 


+  bTu 
I  2 


(4.23) 


The  lower  boundary  condition  Eq .  (4.19)  requires 
W(o)  =  *Jn(o  )kkl 

=  -  Un(o)  Hj  +  X.Un(o)HR 

=  WR(o)  +  <Wj(o) 

where  WR(o)  =  aRUj(o)  +  bRu2 (o)  =  -  Un(o)HI 
and  Wj(o)  =  ajU^o)  +  bjU2(o)  =  Un(o)HR  . 

Case  A:  (k2  -  L2^  z  >  0  . 


(4.24) 


In  this  case  the  upper  boundary  condition  is  given  by 
Eq.  (4.20)  where  Eq.  (4.22)  exhibits  an  exponential  solution. 
Then  at  z  =  z^, 


dW 

dz 


y/K2  -  L2‘  W  • 


MR 
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In  terms  of  coefficients  a  and  b,  this  boundary  condition 


translates  to 


bR  =  -  ^  -  L2  aR 


bj  =  -  VK.2  -  L2  aj 


(4.25) 


Pnth  a  and  b  can  be  evaluated  from  Eq .  (4.24)  and  (4.25) 
stituting  in  Eq.  (4.23),  we  then  get  the  solution 


WR  =  -  u(z)Uli(o)HI/u(o) 


Wj  =  u(z)Un(o)HR/u (o) 


(4.26) 


where  u(z)  =  u^zJ-lK2  -  L2 


'z-z?  n2(Z) 


The  quantity  u(z)  is  obtained  by  solving  the  following  initial 
value  problem: 

Gzz  +  (K2  -  l2)  5  -  0 


wit  h 


“(zt)  ■ 


an“  “z(zt)  ■-  (KJ  -  lTz-zt  ■ 

Case  B:  (K2  -  Lz|  <  0  . 

\  /  ZT 

For  this  case  the  upper  boundary  condition  is  given  by 

£  - z  ss"(us)(L!  -  k2)*  ■ 


(4.27) 
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In  terms  of  coefficients  a  and  b,  this  condition  can  be  written 


ac 


bR  =  -  Sen(us)A2  -  a, 
bj  =  Sgn(Us)>/L2  -  K2  aR  . 


(4.28) 


As  before,  a  and  b  can  be  evaluated  from  Eq.  (4.24)  and  (4.28) 
and  then  substituted  in  Eq.  (4.23)  to  give  the  following  solution 


| 

WR  =  Xoj 


XtU2(o)Hr  -  u i ( o)Hj 


Uj(z) 


^TUi ( ° )HR  +  ^TU  2 ( ° ) H ! 


u  (2)! 


i 

W!  =  Xo  [' 


ui (°)Hr  +  Xtu2(o)Hi 


Uj(z) 


Xtu2(o)Hr  -  Xtu1(o)Hi 


u0(z)| 


Where  X  = 


U. 


o  u2  +  (L2  -  K2 )u2 

xt  -  sAvMvi  - 


2  _  k2)  ^ 


z  =  0 


z  =  zr 


(4.29) 


Unlike  Case  A,  here  we  have  to  solve  two  initial  value  problems 
for  ux  and  u2  subject  to  initial  values  given  by  Eq.  (4.22).  In 
actual  computation  z^,  is  taken  at  the  top  of  the  computational 
region  unless  a  critical  level  (where  Us  vanishes)  is  detected 
at  a  lower  altitude,  in  which  case  zT  is  taken  at  the  closest 
vertical  mesh  below  the  critical  layer. and  reflection  or  absorp¬ 
tion  condition  is  applied  there.  Integration  of  Eq .  (4.16)  then 
ranges  only  between  the  ground  and  the  critical  level. 


4.2  LINMET  SIMULATIONS 

As  a  first  step  LINMET  calculations  are  compared  to 
simple  analytic  solutions  of  the  differential  equations 
Eq .  (4.14).  Under  neutral  atmospheric  conditions  with  constant 
horizontal  velocity  U  ,  the  solution  to  Eq .  (4.14)  in  two- 
dimansions  subject  to  boundary  conditions  similar  to  those  given 
bv  Eq .  (4.19)  and  (4.20)  can  be  obtained  as: 


W(z)  =  -  tk^UoH  exp 


(4.30) 


and  the  horizontal  perturbation  velocity  can  be  obtained  as 


u(z)  =  -  kxUQH  exp  (- kxzj 


(4.31) 


Table  4.1  shows  the  computed  vs.  the  analytic  results  for  the 
horizontal  surface  velocity  u(o)  for  different  modes.  The  com¬ 
parison  is  good, and  the  difference  is  mainly  because  of  discreti¬ 
zation.  Note  that  the  analytic  solution  is  for  the  differential 
equation,  Eq .  (4.14),  whereas  the  LINMET  solution  is  based  on 
the  difference  equation,  Eq .  (4.16). 

Klemp  and  Lilly  (1975)  have  obtained  an  analytic  solution 
of  their  hydrostatic  model  with  constant  stability  and  constant 
velocity  Uo  (one-layer  atmosphere).  They  have  shown  that 


u(o)  = 


(4.32) 


■  (-  f  if  H 

'  o  ' 


LINMET,  which  solves  the  nonhydrostatic  problem,  was  modified  by 
settling  K  ‘  0  to  compare  the  ret  ilts  'with  Eq .  (4.32).  The  com¬ 
puted  vs.  the  analytic  values  are  tabulated  in  Table  4.2.  The 
agreement  is  again  very  close. 


* 


Table  4.1  LINMET  vs.  Analytic  Results  for  u(o). 


Wave  Number 

LINMET 

Result 

Analytic 

Result 

2 

-63.26 

-60.19 

3 

18.22 

17.68 

4 

-  2.08 

-  1.98 

5 

-  0.72 

0.07 

Table  4.2  LINMET  vs.  Analytic  Results  for  u(o) 
for  the  One-Layer  Model . 


Wave  Number 

LINMET 

Result 

Analytic 

Result 

2 

-403.87 

-401.10 

3 

60.46 

60.05 

4 

-  4.67 

-  4.64 

5 

-  1.22 

-  1.22 

As  before,  the  PASS  measurement  records  for  the  day  19 
November  1974  were  used  for  LIHMET  simulations,  since  the  day  was 
characterized  by  relatively  high  winds  and  a  near  neutral  lapse  rate. 
The  structure  of  the  wind  field  in  the  area  is  indicated  by  the  upper 
air  records  of  east-west  and  north-south  wind  components  and  tempera¬ 
ture  presented  for  eight  stations  in  Figures  2.4  through  2.11.  The 
region  of  the  solution  domain  is  indicated  in  Figure  2.15,  which  is 
155  km  wide  in  both  north-south  and  east-west  directions.  As  shown, 
the  area  is  approximated  by  32  x  32  grids,  with  5  km  grid  spacings  in 
both  directions.  The  southwest  corner  is  located  at  local  UTM  coor¬ 
dinates  (340  km,  3570  km).  A  linear  boundary  smoother  is  utilized, 
which  sets  the  boundary  of  the  region  at  the  constant  level  of  1200  m 
MSL  to  avoid  abrupt  discontinuity  from  periodic  boundary  conditions. 
The  top  of  the  computing  mesh  is  set  at  5000  m  MSL,  and  the  bottom  of 
the  computing  region  is  chosen  at  the  mean  height  of  the  terrain, 
which  is  1355  m.  The  vertical  zoning  utilizes  20  levels  with  uni¬ 
form  grid  size  of  192  m.  We  have  not  carried  out  any  sensitivity 
test  with  respect  to  vertical  zoning  and  the  height  of  the  computing 
mesh,  but  we  believe  that  under  certain  conditions  the  solution  will 
be  sensitive  to  finer  zoning  and  a  larger  computing  domain. 

For  the  first  simulation,  LINMET  was  initialized  with  the 
soundings  from  TSX.  The  perturbation  surface  velocity  vector  plot 
is  shown  in  Figure  4.28.  This  is  repeated  for  the  385  m  AGL  and 
3261  m  AGL  in  Figures  4.29  and  4.30,  respectively.  The  contour  of 
normal  velocity  component  (which  is  not  shown  here)  reflects  the 
zero  normal  velocity  boundary  condition,  resulting  in  positive 
vertical  velocity  upslope  on  the  San  Andres  and  Sacramento  Mountains 
and  negative  vertical  velocity  downslope.  At  a  higher  level  (395  m 
above  mean  height  of  the  terrain),  Figure  4.29  shows  higher  distur¬ 
bance  horizontal  velocity  on  the  slopes  but  almost  zero  velocity  in 
the  valleys  where  the  terrain  is  flat.  This  picture  remains  un¬ 
changed  even  at  high  levels  such  as  in  the  ones  shown  in  Figure  4.30, 
which  corresponds  to  3261  m  AGL.  This  is  because  of  the  atmospheric 
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Figure  4.28  Perturbation  Surface  Horizontal  Velocity  Vector 
Plot  for  the  White  Sands  Region.  Flow  Field  was 
Initialized  using  Data  from  Station  TSX  on  11/19/74 
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Figure  4.29  Perturbation  Horizontal  Velocity  Vector  Plot  at 

385. m  above  Mean  Terrain  over  White  Sands  Region. 
Flow  Field  was  Initialized  using  Data  from  Station 
TSX  on  11/19/74. 


Figure  4.30  Perturbation  Horizontal  Velocity  Vector  Plot  at 

3261. m  above  Mean  Terrain  over  White  Sands  Region. 
Flow  Field  was  Initialized  using  Data  from  Station 
TSX  on  11/19/74. 
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conditions  on  that  day,  which  permitted  outgoing  waves  at  the  upper 
boundary.  This  causes  the  disturbance  to  propagate  up  without  any 
reflection.  The  linear  model  obviously  cannot  exhibit  any  low- 
level  turning  that  could  be  expected  because  of  early  morning  low- 
level  inversion  (~200  m)  as  shown  by  the  temperature  data  at  TSX. 

The  second  simulation  was  made  by  initializing  with  measure¬ 
ment  data  from  station  0R0.  The  location  of  0R0  is  close  to  TSX,  and 
both  are  located  in  the  Tularosa  Valley.  Wind  and  temperature  at  0R0 
are  close  to  those  at  TSX,  and  the  simulation  did  not  exhibit  any 
new  trend  or  feature. 

The  third  simulation  was  made  by  initializing  LINMET  with 
data  from  station  RAM.  For  this  simulation,  velocity  vector  plots 
are  shown  for  the  surface  level,  385  m  level  and  3261  m  level  in 
Figures  4.31  through  4.33,  respectively.  The  low  level  flow  in 
the  simulation,  when  superimposed  on  the  ambient  wind,  looks  very 
close  to  the  VARMET  stable  simulation  (presented  in  Chapter  5) 
using  upper  air  data  for  all  eight  stations  (Figure  5.2). 

We  have  encountered  problems  in  initializing  LINMET  with 
WAR,  SMR,  and  HMS  data.  Insufficient  vertical  resolution,  along 
with  possible  occurrence  of  resonant  modes  resulting  in  large 
amplitude  waves  (Klemp  and  Lilly,  1975)  could  have  caused  this. 

The  problem  could  be  resolved  by  employing  finer  vertical  grid 
spacings  and  extending  the  horizontal  domain  of  computation. 

The  limited  area  chosen  from  the  West  German  region  for 
simulation  is  shown  by  Figure  2.14.  The  area  is  approximated  by 
32  x  32  grids  with  1  km  grid  spacing  in  both  directions.  The 
southwest  corner  of  this  area  is  located  at  local  UTM  coordinates 
(526.2  km,  5584.3  km).  The  linear  boundary  smoother  is  applied 
to  a  boundary  strip  5  km  wide,  which  sets  the  boundary  terrain 
heights  at  a  constant  value  of  400  m.  The  top  of  the  computing 
region  is  set  at  5000  m  MSL  and  the  bottom  at  379  m  MSL,  which 
is  the  average  height  of  the  region.  The  vertical  zoning  utilizes 
20  levels  with  uniform  grid  size  of  243  m. 


Figure  4.31  Perturbation  Surface  Horizontal  Velocity  Vector 
Plot  for  the  White  Sands  Region.  Flow  Field  was 
Initialized  using  Data  from  Station  RAM  on  11/19/74 


Figure  4.32  Perturbation  Horizontal  Velocity  Vector  Plot  at 
385  m  above  Mean  Terrain  over  White  Sands  Region. 
Flow  Field  was  Initialised  using  Data  from  Station 
RAM  on  11/19/74. 


Figure  4.33  Perturbation  Horizontal  Velocity  Vector  Plot  at 

3261  m  above  Mean  Terrain  over  White  Sands  Region. 
Flow  Field  was  Initialized  using  Data  from  Station 
RAM  on  11/19/74. 
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Using  the  upper  air  soundings  from  Meiningen,  East  Germany 
(Figure  2.16),  which  is  east  of  the  region  of  simulation,  a  LINMET 
simulation  was  carried  out.  Figure  4.34  shows  the  surface  perturba¬ 
tion  velocity  vector  for  this  simulation.  This  plot  shows  some  oscil¬ 
lations  that  are  believed  to  be  related  to  the  finite  extent  of  the 
horizontal  domain.  Even  though  the  terrain  in  question  is  without 
large  gradients,  sometimes  it  is  necessary  to  extend  the  computing 
domain  to  great  distance  over  flat  terrain  to  effectively  treat  the 
lee  waves,  which  may  otherwise  be  truncated  with  periodic  boundary 
conditions  over  the  limited  terrain.  The  smaller  horizontal  grid 
size  (1000  m)  could  pose  a  limitation  for  hydrostatic  models,  but 
as  LINMET  is  a  nonhydrostatic  model  this  limitation  does  not  apply. 
More  experiments  and  sensitivity  studies  must  be  carried  out  before 
definite  conclusions  can  be  reached. 
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Figure  4.34  Perturbation  Surface  Horizontal  Velocity  Vector 

Plot  for  the  West  Germany  Region  containing  Fulda. 
Flow  Field  was  Initialized  using  Upper  Air  Data 
from  Meiningen,  East  Germany. 


5.  OBJECTIVE  ANALYSIS  MODEL  VARMET 


In  this  category  of  simplified  models,  one  follows  a  pro¬ 
cedure  of  interpolation  and  extrapolation  first  to  determine  an 
estimated  wind  field  from  observed  data.  This  field  is  then 
adjusted  to  account  for  terrain  effects  and  atmospheric  stability 
considerations  constrained  by  the  condition  that  the  resulting 
wind  field  be  nondi vergent .  In  addition  to  the  mass  conservation 
constraint,  one  can  impose  momentum  and  energy  conservation 
constraints  as  well,  but  at  the  cost  of  added  complication  and 
computing  expense.  The  accuracy  so  gained  may  be  lost  if 
there  is  considerable  uncertainty  in  the  observed  data.  In  the 
formulation  to  follow,  we  have  therefore  restricted  our  work  to 
the  mass  conserving  wind-field  models. 

5.1  VARMET  FORMULATIONS 

The  theoretical  basis  for  the  model  was  developed  by 
Sasaki  (1970)  and  later  used  in  this  context  by  Sherman  (1977). 
For  this  model  the  squared  variation  of  the  wind  field  is  mini¬ 
mized  subject  to  the  constraint  that  the  adjusted  field  be  non- 
divergent.  The  irregular  surface  boundary  z  =  zg(x,y)  is  first 
removed  for  computational  purposes  by  a  coordinate  transformation 
in  which  the  surface  becomes  a  coordinate  surface  as  well.  It 
is  convenient  to  c^rform  a  nonorthogonal  transformation  in  which 
only  the  vertical  coordinate  z  is  replaced  by  a  function  a  , 
which  depends  on  x ,  y  as  well  as  z.  The  new  coordinate  is  taken 
linear  in  z: 

Zfp— z  z»-z 

_  _  i  I  /  R  1  \ 


where  z T  is  the  constant  altitude  of  the  top  mesh  and  zg  is 
the  terrain  surface.  The  values  of  o'  are  in  the  range  of 
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0£ o  il  (a= 0  corresponds  to  the  top  of  the  mesh,  o=l  corresponds 
to  the  surface).  The  variation  integral  written  in  the  x-y-z 
frame  is 


[(u"u° )  +(v-ve)]+  a£  (w-w0 ) 


(5.2) 


where  c^,  a?  are  the  gauss  precision  moduli  (proportional 
to  1/(2B2)  ,  B  being  observation  error  and/or  deviation  of 
the  observed  field  from  the  desired  adjusted  field),  and  X  is 
the  Lagrange  multiplier.  The  observed  field  is  represented  by 
u0  ,  v0  ,  and  w0  .  The  weights  a  may  be  parametrized,  that 
is,  constants  for  a  given  trial  but  may  vary  from  trial  to  trial. 
More  generally,  the  weights  should  be  a  function  of  position 
within  the  region,  being  large  in  the  vicinity  of  the  actual 
observation  and  small  in  regions  remote  from  observations.  The 
ratio  (  ai/a2  )  can  be  viewed  as  a  stability  parameter  also. 

Since  in  a  strongly  stable  flow,  the  vertical  motion  is  inhibited, 
(  ai/a2  )  should  be  assigned  a  low  value  for  stable  ambient  con¬ 
dition.  The  ratio  equal  to  one  should  be  chosen  for  neutral 
flow  condition  since  there  is  no  preferential  direction  for  motion 
under  those  conditions.  The  associated  Euler-Lagrange  equations 
to  minimize  I  can  be  written  as 


2ct12(u-u0) 
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From  the  nature  of  the  coordinate  transformation  in  Eq .  (5.1) 
we  can  derive  expressions  for  partial  derivatives  of  an  arbi¬ 
trary  function  f  in  the  two  representations. 
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The  subscripts  on  partial  derivatives  indicate  which 
variables  are  being  held  constant  during  differentiation,  and 
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Using  this  coordinate  transformation,  the  Euler-Lagrange  equation, 
Eq .  (5.3),  can  be  rearranged  as  follows: 


=  v0  + 


w  =  w0  - 
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(5.5) 


Also  using  the  transformation/  the  divergence-free  condition 

V-u  =  0  ,  or 


3_u  ,  _3v  ,  3w  n 
3x  3y  3z  "  u 
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takes  the  form  in  o  -  coordinate 


3ttu  3ttv  3ttw  _ 

3x  3y  3o 


(5.6) 


where , 

ttw  =  a  (zsxu  +  zSyv  )  -  w 


(5.7) 


The  Eq.  (5.6)  is  written  in  the  so-called  "conservative  form," 
in  which  flux  analogy  can  be  directly  applied.  This  form  is 
useful  in  forming  difference  equations  which  are  more  accurate 
through  being  constrained  to  conserve  mass  in  the  difference 
equations . 

We  now  can  eliminate  u,  v,  and  w  in  Eq.  (5.6)  by  using 
Eq.  (5.5)  and  (5.7)  to  arrive  at  the  following  equation  in  X  . 
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The  boundary  conditions  for  the  above  consistent  with  minimization 
of  I  are  given  by 

A6(un)  =  0  (5.9) 

where  6(un)  is  variation  of  the  velocity  component  normal 
to  the  boundary.  When  X  is  zero  on  the  boundary,  its  normal 
derivative,  in  general,  is  not  zero  and  results  in  a  new  normal 
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boundary  velocity  different  from  the  initial  guessed  value. 
Therefore  X  =  0  is  appropriate  for  "open  boundaries."  On  the 
other  hand,  when  un  is  held  fixed,  the  Euler-Lagrange  equations 
imply  that  the  normal  derivative  of  X  should  vanish.  This  may  be 
appropriate  for  "no-f low-through"  boundaries.  Clearly,  the  latter 
is  applicable  to  top  and  bottom  boundaries  where  w  =  0  .  On 
the  lateral  sides,  either  kind  of  boundary  conditions  may  apply, 
largely  depending  on  the  density  of  observation  data  near  the 
boundaries.  Observation  data  are  usually  available  at  a  few 
stations.  These  values  should  be  judiciously  used  to  estimate 
u0  ,  vo, and  w0  at  all  grid  points  through  a  process  of  inter¬ 
polation/extrapolation.  One  such  method,  which  is  presently 
being  used,  is  described  in  Section  2. 

At  present  on  all  four  lateral  sides 


X  =  o  (5.10) 

is  used  At  top  and  bottom  w  =  0  or  w  =  c(zSxu  +  Zg^v)  con¬ 

dition  is  used.  Using  Eq  .  (5.5)  and  (5.7),  this  translates  to 
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at  the  bottom. 

The  Eq .  (5.8),  which  is  a  Helmholtz's  equation  for  X 
along  with  the  boundary  conditions  Eq.  (5.10)  and  (5.11),  com¬ 
pletes  formulation  of  the  problem.  In  practice,  iterative  methods 
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making  use  of  over-relaxation  technique  are  employed  to  solve 
the  system.  The  numerical  procedure  is  described  below. 


A  staggered  grid  network  is  used  to  discretize  the 
Helmholtz  Equation,  Eq.  (5.8),  for  numerical  solution.  In  a 
cube  with  sides  Ax,  Ay,  and  Ac  ,  the  u,  v,  and  w  components  of 
the  velocity  are  defined  at  the  center  of  the  faces  normal  to  x, 
y  ,  and  o  directions,  respectively.  The  Lagrange  multiplier 
A  is  defined  at  the  geometric  center  of  the  cube.  With  this 
staggering,  one  can  write  a  very  compact  finite  difference  equa¬ 
tion  approximating  Eq.  (5.8)  with  minimum  truncation.  The 
equation  can  be  solved  by  the  method  of  Successive-Point-Over- 
Relaxation  with  an  optimum  value  for  successive  parameter  w. 

The  rate  of  convergence  strongly  depends  on  the  aspect  ratio 
of  the  mesh,  modified  by  04  and  a 2  ,  namely, 

a  l2 

Max  (Ax, Ay)  1 

tt  .  Aa  .  a„ 
min  mm  2J 

Typically,  this  value  is  very  large  compared  to  unity  (partly 
because  vertical  zone  size  is  much  smaller  than  the  horizontal 
zone),  resulting  in  no  or  very  slow  convergence.  Other  investi¬ 
gators,  while  following  such  formulation,  have  been  plagued  by  the 
same  problems  in  the  past.  We,  therefore,  made  use  of  Line-Over- 
Relaxation  in  which  the  A  values  in  a  given  vertical  column  are 
calculated  simultaneously  rather  than  successively.  Because  only 
terms  involving  second  or  lower  derivatives  with  respect  to  a 
appear  in  Eq.  (5.8),  the  implicit  treatment  in  a-direction  results 
a  system  of  tridiagonal  equations  in  the  discretized  space,  one 
for  each  column.  These  can  be  readily  solved  using  the  well-known 
Gaussian  elmination  procedure.  This  method  of  solution  for  Eq. 
(5.8)  converges  reasonably  fast  to  yield  a  divergence-free  solu¬ 
tion  for  the  velocity  field,  and  addin'g  very  little  to  overall 
computation  expense. 


The  implicit  line  relaxation  method  in  VARMET,  while  making 
the  solution  insensitive  convergence-wise  to  physical  modeling 
parameters  such  as  a } / u 2  and  surface  slopes,  shows  instability 
for  small  values  of  ^\la2  (  '  1  / 01  >  <  0-01)-  Since  c^/a^,  -*■  0 
decouples  the  equations  in  the  vertical  direction,  a  situation 
typical  under  strong  stable  atmospheric  condition,  the  resulting 
reduction  in  the  dimensionality  of  the  problem  can  possibly  explain 
this  instability  in  the  relaxation  scheme.  In  any  event,  we  have 
found  that  the  addition  of  an  explicit  damping  term  stabilizes  the 
solution  in  such  cases.  The  form  of  the  damping, 


was  chosen  to  be  added  to  Eq.  (5.8)  based  on  the  magnitude  of  the 
terms  in  the  difference  expressions  used  in  the  iterative  scheme. 

It  has  been  found  that  a  value  of  f  ^  0.05  to  0.10  provides  suffi¬ 
cient  stabilization  while  retaining  acceptable  convergence.  This 
may  vary  from  problem  to  problem;  however,  no  general  parameteri¬ 
zation  has  been  developed  to  date. 

5.2  VARMET  SIMULATIONS 

The  procedure  is  applied  to  the  wind  data  obtained  from 
PASS  experiment  on  19  November  over  the  limited  White  Sands  terrain 
def  i  .ied  by  the  dashed  lines  in  Figure  3.3.  The  coefficients  a! 
and  a-,  ,  which  are  believed  to  be  stability-dependent,  but  the 
precise  nature  of  which  is  yet  to  be  evaluated,  are  set  equal  to 
each  other.  The  over-relaxation  parameter  oj  is  set  to  a  constant 
value  of  1.7.  In  the  east-west  directions  38  grid  points,  and  in 
the  north-south  directions  32  grid  points  ,  are  employed  with  an 
equal  grid  interval  of  25  km.  The  southwest  corner  (as  shown 
by  the  dashed  lines  in  Figure  3.3)  is  .located  at  the  local  UTM 
coordinates,  337.5  km  east  and  3562.5  km  north.  In  the  vertical, 
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12  zones  are  used  with  nonuniform  zone  thickness  defined  by  a 
geometric  progression  with  1.37  expansion  ratio,  the  one  nearest 
to  the  surface  being  approximately  30  m  and  the  farthest  one 
being  1000  m.  The  top  of  the  computing  mesh  is  set  at  a  constant 
altitude  of  4800  m  MSL, and  the  convergence  criterion  for  >.  is 
set  at  10~-  for  all  the  simulations.  The  program  makes  use  of 
a  graphics  package  to  prepare  vector  and  contour  plots  of  various 
quantities  of  interest  at  user-specified  intervals  on  microfiche. 

A  typical  run  on  this  grid  takes  about  80-100  iterations  to 
converge  and  uses  about  40  CP  sec  on  a  CDC  7600  computer,  which 
includes  the  time  spent  in  graphics. 

In  an  effort  to  carry  out  a  systematic  verification  over 
the  White  Sands  region,  we  have  built  a  simulation  matrix  designed 
to  investigate  the  efficacy  of  the  extrapolation  technique  for 
utilizing  surface  station  wind  data  as  well  as  the  overall  simula¬ 
tion  accuracy  with  various  numbers  of  data  stations.  Figure  5.1 
defines  the  10  cases  that  make  up  the  matrix.  In  the  matrix,  u 
indicates  that  the  actual  air  record  (interpolated  to  the  finite 
difference  grid  locations)  was  used  in  the  wind-field  initialization, 
s  indicates  that  the  data  were  used  as  a  surface  station  (see  below), 
and  x  indicates  that  the  data  were  unused  in  the  initialization  and 
is  thus  available  for  comparison  with  the  calculated  results. 

The  use  of  surface  stations  requires  some  explanation, 
since  in  most  applications  of  the  model  only  near-surface  winds 
will  be  available.  In  the  usual  treatment  of  such  data  in 
VARMET,  the  near-surface  wind  is  extrapolated  in  the  vertical 
as  follows.  The  wind  is  assumed  to  take  on  an  equilibrium 
boundary  layer  profile  with  power-law  6  (which  may  depend  on 
terrain,  surface  roughness,  stability,  etc.)  up  to  an  arbitrary 
boundary  layer  height  [0(200  m)],  then  a  linear  interpolation  of 
wind  speed  and  angle  is  used  up  to  some  assumed  uniform  high  level 
winds  at  the  top  of  the  computational  *mesh.  In  the  present  calcu¬ 
lations,  this  procedure  was  used  by  taking  the  lowest  value  in  the 
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Data  Station 


Figure  5.1  Matrix  of  Verification  Simulations 
for  11/19/74, 
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upper  air  record  (usually  ^200  m  above  ground)  and  extrapolating 
within  a  200  ra  AGL  boundary  layer  with  6  =  1/7  and  extrapolated 
above  the  boundary  layer  using  an  assumed  average  value  for  the 
high  level  wind.  It  should  be  mentioned  that  the  200  m  value  for 
the  boundary  layer  thickness  is  typical  of  what  would  be  assumed 
if  nothing  were  known  about  the  wind  profile  so  that  it  does  not 
presuppose  knowledge  of  the  upper  air  records.  Also,  the  values 
of  the  upper  level  wind  speed/direction  were  taken  as  a  rough 
average  of  all  the  upper  air  soundings,  but  an  equally  representative 
value  could  be  determined  from  the  gradient  wind  approximation  to 
synoptic  weather  data  in  the  general  case.  Thus,  it  is  not  believed 
that  this  choice  will  bias  the  results. 

Cases  2,  4,  and  6  in  the  test  matrix  were  chosen  to  compare 
overall  flow  results  utilizing  4,  6  and  7  stations  taken  as  surface 
stations,  respectively.  Case  1,  of  course,  used  all  of  the  upper 
air  records  and  represents  the  "best  estimate"  of  the  wind  field 
for  the  19  November  1974  event  and  therefore  serves  as  a  basis  for 
comparison  for  the  other  cases.  Wind  field  vector  and  speed  contour 
plots  are  presented  for  Cases  1,  2,  4,  and  6  in  Figures  5.2 
through  5.5,  respectively,  for  the  near  surface  (15  m  AGL)  flow. 

The  heavy  arrows  in  each  vector  plot  correspond  to  the  original  data 
at  each  station  interpolated  to  the  15  m  AGL  level.  Comparison  of 
the  results  shows  very  good  qualitative  and  quantitative  correspon¬ 
dence  of  the  flow  in  each  case.  Even  Case  6,  which  uses  only  one 
actual  upper  air  record  (station  SMR)  and  seven  surface  station 
prescriptions,  compares  well  with  the  Case  1  baseline  result, 
especially  in  terms  of  magnitude  of  wind  speed,  throughout  the  solu¬ 
tion  region. 

It  should  be  reiterated  that  the  surface  station  prescrip¬ 
tion  used  the  single  element  of  the  raw  data  record  that  corresponded 
to  ^200  m  AGL.  In  the  data  record  for  RAM,  there  was  substantial 
variation  in  the  data  below  this  level.’  Thus,  it  is  not  surprising 
that  the  simulations  using  surface  stations  miss  some  details  of 


Figure  5.2  Near-Surface  Wind  Field  Velocity  Vector 
and  Speed  Contour  Plots  (15  m  AGL),  Case 
1  (all  stations  used  as  upper  air  records) 
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Figure  5.3  Near-Surface  Wind  Field  Velocity  Vector 
and  Speed  Contour  Plots  (15  m  AGL) ,  Case 
2  (stations  TSX  (1),  WAR  (3),  RAM  (6), 
and  APA  (7)  used  as  surface  stations  in 
initialization) . 
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Figure  5.4 


Near  Surface  Wind  Field  Velocity  Vector 
and  Speed  Contour  Plots  (15  m  AGL),  Case 
4  (all  stations  except  SMR  (4)  and  HMS 
(8)  used  as  surface  stations  in  initial¬ 
ization)  . 
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Figure  5.5  Near  Surface  Wind  'Field  Velocity  Vector 

and  Speed  Contour  Plots  (15  m  AGL),  Case 
6  (all  stations  except  SMR  (4)  used  as 
surface  stations  in  initialization). 


the  near-surface  flow  shown  in  the  Case  1  result,  for  example, 
the  eddy-like  structure  in  the  valley  encompassing  the  nearby 
stations  0R0  (westerly  flow)  and  RAM  (easterly  flow)  and  the 
reverse  flow  indicated  at  station  WAR.  These  items  notwithstand¬ 
ing,  the  simulated  results  with  the  surface  station  prescription, 
used  in  lieu  of  the  detailed  upper  air  records,  compare  very  well 
with  the  baseline  flow  that  utilizes  all  the  upper  air  records. 

The  other  cases  in  the  matrix  can  be  used  to  verify  the 
accuracy  of  the  model  predictions  at  unused  data  stations.  Results 
of  wind  speed  scatter  plots  (q„_-._  vs.  q  ,  )  are  thus  presented  for 

Ca  J.  C  O DS 

Cases  1,  3,  5,  9,  and  10  in  Figures  5.6  through  5.10,  respectively. 

In  each  figure,  the  results  at  15  m,  250  m,  500  m,  and  1000  m  AGL 
are  compared  with  the  corresponding  upper  air  record  data.  In  Case 
1  (Figure  5.6),  of  course,  all  the  data  were  used  in  the  initial¬ 
ization  process,  so  it  is  not  surprising  that  the  calculated  wind 
speeds  compare  so  well  with  the  observed  wind  speeds.  It  is  pre¬ 
sented,  however,  to  show  that  there  is  little  "drift"  from  the 
initialized  values,  at  least  in  regions  of  flat  terrain,  which 
characterizes  all  of  the  data  stations.  It  also  indicates  the 
optimum  degree  of  comparison  that  can  be  expected  for  the  unused 
data  stations. 

The  scatter  plot  of  the  unused  data  stations  indicate 
substantial  scatter  in  .  o  predictions  and  show  that  Case  5  (Figure 
5.8),  which  utilizes  only  two  upper  air  records  in  the  initializa¬ 
tion,  compares  as  well  to  the  observations  as  the  other  cases  that 
utilize  more  initialization  stations.  In  each  case,  the  simulations 
show  a  clear  tendency  to  overpredict  the  wind  speeds  in  the  valley 
(where  all  the  data  stations  occur).  However,  the  winds  in  the 
high  speed  areas  (mountainous  regions)  are  in  relatively  good  cor¬ 
respondence  for  each  case.  This  tendency  to  overpredict  the  winds 
can  be  traced,  at  least  in  part,  to  the  uniquely  high  (not  present 
in  the  other  records)  and  low  level  wirrds  at  station  SMR.  This  is  sub 
stantiated  by  the  Case  10  simulation,  which  is  the  only  case  that 
does  not  use  SMR  in  the  initialization.  As  shown  in  Figure  5.10, 
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Figure  5.6  Comparison  of  Calculated  to  Observed  Wind  Speeds: 

WSMR,  11/19/74,  Case  1  (all  data  stations  used 
in  initialization). 
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the  predicted  results  for  Case  10  are  quite  good  (lower  speeds) 
for  all  stations  except  SMR.  Unfortunately,  for  this  choice  of 
a  comparison  matrix,  all  other  cases  used  SMR  in  the  initialization. 

The  comparisons  presented  here  provide  clear  verification 
of  the  initialization  procedure  for  surface  station  data  indicating 
that  the  effect  of  the  vertical  structure  at  the  initialization 
stations  is  not  of  overriding  importance, especially  in  relation  to 
the  prediction  of  the  magnitude  and  extent  of  high  wind  speed  areas. 
Moreover,  they  have  shed  some  light  on  the  accuracy  of  the  results 
as  influenced  by  the  number  of  initialization  stations  and  indicate 
that  anomalous  data  (station  SMR)  can  have  a  biasing  effect  on  the 
simulation  results. 


6.  SUMMARY  AND  CONCLUSION 


The  objective  of  this  study  is  to  identify,  develop, and 
test  improved  simulation  techniques  for  mesoscale  meteorology. 

The  emphasis  is  placed  on  flow  over  terrain  features  with  the  goal 
of  obtaining  one  or  more  computer  models  that  are  at  the  same 
time  reasonably  accurate,  easy  to  use,  and  economical  of  computing 
expense.  This  has  been  achieved  in  the  study;  in  fact,  one  code 
was  adapted  for  this  study,  and  two  other  codes  were  developed  from 
fundamentals . 

The  primitive  equation  SIGMET  code  was  adapted  from  a  pre¬ 
vious  version.  This  model,  with  an  advanced  prescription  of  boundary 
layer  turbulence,  radiation,  soil/air  interaction,  and  moisture 
treatment  was  used  to  simulate  atmospheric  events  over  the  White 
Sands  Region  for  selected  days  in  the  period  November  -  December  1974 
The  results  compared  favorably  with  the  PASS  data  obtained  during 
the  same  period  in  the  region. 

An  extensive  survey  which  included  a  review  of  recent  lit¬ 
erature  on  orographic  flow  and  discussion,  with  several  experts  in 
mesoscale  modeling,  was  undertaken  in  Phase  I  of  the  study  in  order 
to  recommend  a  simplified  physics  model  for  further  development  and 
testing.  This  survey  has  been  reported  in  detail  in  the  Phase  I 
report  (Patnaik,  Freeman;  1977).  The  two  models  that  emerged  as 
a  result  of  this  survey  are  VARMET  and  LINMET.  The  formulation  of 
VARMET  is  based  on  the  objective  analysis  approach  to  wind  field 
synthesis.  The  model  uses  terrain  conformal  coordinate  for  accurate 
representation  of  terrain  effects.  The  degree  of  success,  using  the 
code,  depends  heavily  on  the  number  and  location  of  the  data  stations 
that  are  used  to  initialize  the  calculation.  For  data-rich  regions 
the  VARMET  simulations  are  very  close  to  real  events, whereas, for 
data-sparse  regions, the  simulations  ard  rather  poor. 

A  formulation  containing  more  fundamental  principles  than 
VARMET, but  less  sophisticated  than  SIGMET,  resulted  in  the  code  LINMET 
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It  is  based  on  the  3-D,  steady-state,  linearized,  inviscid  Navier- 
Stokes  equations.  Calculations  using  LINMET  are  not  very  expensive, 
but  it  does  simulate  realistic  flow  features  based  on  very  little 
data.  Unfortunately,  only  limited  experiments  have  been  carried 
out  using  LINMET,  which  are  reported  in  this  study. 

A  considerable  amount  was  learned  about  orographic  flow  and 
mesoscale  modeling  in  this  study.  Advantages  and  shortcomings  of 
various  classes  of  models  were  identified  and  understood.  The 
linear  code  LINMET,  in  particular,  showed  very  promising  features 
which  warrant  further  investigation  in  the  future.  Addition  of  a 
boundary  layer  treatment  and  a  terrain  conformal  coordinate  system 
(as  in  SIGMET  and  VARMET)  to  LINMET  will  improve  its  performance 
considerably  while  keeping  the  computational  effort  to  a  minimum; 
this  would  then  prove  to  be  a  viable  alternative  to  SIGMET.  It  is 
also  possible  to  include  certain  aspects  of  nonlinearity  to  LINMET 
following  the  work  of  Smith,  1977. 
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APPENDIX  A 
SIGMET  ORGANIZATION 


A . 1  SIGMET  ORGANIZATION 

SIGMET  is  written  in  modular  form,  through  the  use  of  sub¬ 
routines,  with  each  module  representing  an  aspect  of  the  physical 
model.  There  is  thus  a  close  correspondence  between  the  sub¬ 
routine  structure  and  the  individual  physical  processes  being  sim¬ 
ulated.  The  logical  flow  of  the  program  is  controlled  by  subprogram 
MAIN  and  proceeds  through  the  subroutines  as  indicated  in  the  summary 
flowchart  of  Figure  A.l. 

The  numerical  time  and  3-D  space  dependent  solution  of  the 
system  of  equations  proceed  sequentially  through  the  subroutines 
indicated  in  the  flow  diagram.  The  solution  state  is  advanced 
forward  in  time  in  a  step-by-step  fashion  by  the  execution  of  the 
subroutines  and  thereby  mimicks  the  actual  evolution  of  the  atmo¬ 
spheric  properties.  One  computational  cycle  advances  all  of  the 
state  variables  by  a  small  increment  corresponding  to  the  intensity 
of  the  physical  processes  calculated  at  the  time  in  question.  The 
increment  in  time  represented  by  one  cycle  is  quite  small  so  that 
the  rates  of  change  of  the  atmospheric  variables  may  be  considered 
constant  over  the  interval.  Thus,  beginning  with  an  initial  con¬ 
dition  state,  the  desired  solution  evolves  iteratively  in  response 
to  physical  processes  and  applied  boundary  conditions. 

Within  each  cycle,  a  sequence  of  computational  steps  is  per¬ 
formed  to  advance  the  solution  through  one  time  interval.  These 
steps  are  performed  in  the  subroutines  indicated  in  Figure  A.l  and 
in  the  order  shown.  Because  of  the  solution's  explicit  nature, 
the  order  of  execution  of  some  of  these  routines  is  not  important; 
the  order  has  been  arranged  primarily  for  programming  simplicity 
and  for  reducing  the  amount  of  computei*  storage  required  to  contain 
the  current  values  of  the  state  variables.  Because  of  the  sequential 
nature  of  the  solution  it  is  not  necessary  to  store  these  variables 
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Figure  A.l  SIGMET  Logical  Flow  Diagram  (all  subroutines 
controlled  and  called  by  MAIN). 


as  a  function  of  time.  Only  instantaneous  values  are  required  to 
continue  the  solution, and  information  regarding  the  solution  state 
at  any  time  is  provided  by  periodic  edits  of  these  values. 


A. 2  SUBROUTINE  DESCRIPTION 

The  subroutines  that  affect  the  numerical  solution  are  now 
described  in  some  detail.  This  includes  a  description  of  the  func¬ 
tion  and  operation  of  each  subroutine  as  well  as  some  additional 
explanatory  comments  concerning  the  physical  model  and  finite  '■'■’t- 
ference  technique.  A  spatially  staggered  mesh  is  utilized  in  e 
model  (Figure  A. 2)  that  is  referenced  in  the  description  to  ;  ’ow. 


X  *  u,  v,  itu,  wv  points 
*  2 

A  *  a ,  q  points 
□  ■  7T,  T,  C,  <p  points 


Figure  A. 2  A  Three-Dimensional  Grid  Element. 
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The  INITAL  subroutine  processes  data,  as  provided  by  data 
cards,  so  that  all  information  required  for  subsequent  normal 
execution  of  a  computation  cycle  is  available.  The  data  cards 
provide  values  for  the  following  quantities: 

•  spatial  increments  of  the  horizontal  and  vertical 
distance 

•  time  intervals 

•  number  of  cycles  of  calculation 

•  terrain-height  array 

•  latitude  and  longitude  of  the  region 

•  number  of  parameters  characterizing  the  terrain  such 
as  volumetric  fractions  of  clay,  sand,  and  organic 
material  in  the  terrain,  roughness  height,  vegetation 
fraction,  etc. 

•  a  number  of  parameters  characterizing  the  solar  and 
terrestrial  radiation  through  the  atmosphere 

•  tables  in  pressure  coordinates  describing  the  vertical 
profiles  of  the  ambient  values  of  horizontal  velocity 
components,  temperature  and  amount  of  water  vapor  in 
the  atmosphere 

•  initial  surface  temperature 

•  initial  surface  wetness  and 

•  indicators  controlling  the  output  of  information 

Options  are  provided  so  that  some  of  the  read-in  quantities  can  be 
overwritten  and  suitable  data  generated  internally  in  INITAL.  For 
example,  the  vertical  coordinate  a  can  be  generated  using  a  log- 
linear  formula  or  a  geometric  progression  series  resulting  in  smooth 
nonuniform  distribution  of  o  levels.  Another  example  is  the 
ambient  temperature  table  can  be  constructed  using  a  temperature 
near  the  surface,  two  lapse  rates, and  an  inversion  pressure  height. 
The  surface  temperature  can  be  defined  using  the  temperature  table, 
and  the  surface  pressure  and  the  water  vapor  concentration  can  be 
calculated  from  relative  humidity  values  instead  of  being  read-in 
directly . 


The  initial  data  are  processed  so  that  they  are  acceptable 
to  the  other  subroutines  of  SIGMET.  In  some  cases,  this  involves 
conversion  of  units  from  those  in  which  the  input  data  are  commonly 
available  to  those  consistent  with  the  formulation.  In  addition,  a 
number  of  auxiliary  quantities  are  formed  that  are  used  during  the 
calculation  but  do  not  change  with  time.  Precalculation  of  these 
quantities  saves  computing  time.  In  other  cases,  the  data  are 
interpolated  to  provide  values  at  locations  specified  by  the  compu¬ 
tational  grid.  This  second  function  is  performed  in  a  subroutine 
called  NESTAR,  which  forms  such  quantities  as  the  pressure  thickness 
(of  the  part  of  the  atmosphere  for  which  calculation  is  performed), 
the  horizontal  velocity  components,  the  temperature,  and  the  initial 
water  vapor  concentration.  These  quantities  are  assigned  to  the 
appropriate  storage  arrays. 

NESTAR 

In  NESTAR  a  pressure  surface  is  calculated  corresponding  to 
a  prespecified  constant  boundary-layer  height  above  the  local  ter¬ 
rain.  Above  this  surface,  the  values  of  the  horizontal  velocity 
components  are  interpolated  from  the  tables;  below  the  surface, 
the  values  are  computed  from  a  power-law  profile.  Temperature  and 
moisture  are  interpolated  from  this  table  at  all  levels.  All  the 
interpolation  is  done  by  calling  the  subroutine  function  ENTERP. 
Subroutine  NESTAR  also  defines  initial  values  of  the  turbulent 
kinetic-energy  arrays  using  the  properties  of  the  ambient  atmos¬ 
phere  and  the  "level  2"  formulation  of  Yamada-Mel lor  turbulence 
prescription . 

FRICH 

The  FRICH  subroutine  performs  two  functions:  horizontal  dif¬ 
fusion  coefficients  are  calculated  as  functions  of  the  local  state 
variables,  and  (using  these  calculations)  the  turbulent  horizontal 
diffusion  taking  place  in  one  time  interval  is  calculated.  The 
result  of  these  steps  is  to  evaluate  the  contribution  of  horizontal 
diffusion  to  the  rates  of  change  of  all  of  the  primary  variables. 
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The  horizontal  diffusion  coefficients  are  obtained  from 
Smith  and  Howard  (1972)  and  are  incorporated  into  FRICH  in  the  form 
of  tables.  Entries  of  the  main  table  are  values  of  the  diffusion 
coefficient  as  a  function  of  altitude  for  each  of  the  standard 
stability  categories  according  to  the  AEC  classification.  A  factor 
accounting  for  the  dependence  on  stability  of  the  ratio  of  horizon- 
tal-to-vertical-dif fusion  coefficient  is  also  applied. 

Because  the  horizontal  diffusion  is  a  slow  process,  the  cal¬ 
culation  can  be  performed  explicitly  without  regard  for  the  time 
interval.  Consequently,  the  second  task  of  the  FRICH  subroutine  is 
to  perform  a  calculation  of  the  increment  to  the  primary  variables, 
u  ,  v  ,  o,and  C  ,  resulting  from  horizontal  diffusion.  Numerical 
stability  of  this  diffusion  calculation  applied  to  the  leapfrog 
method  requires  that  a  special  treatment  of  the  time  centering  be 
used.  If  this  term  were  the  only  one  being  calculated,  a  typical 
equation  would  be  evaluated  as  follows: 
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which  has  the  stability  requirement 


(  A .  1 ) 


(A. 2) 


The  time  centering  of  the  right-hand  term  is  not  permitted,  because 
such  a  scheme  is  numerically  unstable.  Because  of  the  large  hori¬ 
zontal-space  increments,  the  stability  condition  given  by  Eq.  (A. 2) 
is  considerably  less  stringent  than  that  for  gravity  waves  to  be 
discussed  below. 
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FIELDS 


The  FIELDS  subroutine  is  executed  after  FRICH  at  the  begin¬ 
ning  of  each  cycle  to  calculate  auxiliary  quantities,  including 
the  two  components  of  horizontal  velocity,  temperature,  and  moisture. 
These  quantities  are  obtained  from  the  primary  variables,  iru  ,  nv  , 
TtO  ,  and  nC  ,  where  n  *  Pg  -  is  the  pressure  thickness  of  the 
computational  region,  and  0  is  the  potential  temperature.  By 
placing  this  subroutine  after  FRICH,  storage  of  u  ,  v  ,  0,  and  C 
at  two  time  levels  is  avoided.  FRICH  uses  the  old  values  of  these 
variables  for  stability  and  then  FIELDS  updates  them  to  current  time. 

TENDH 

In  TENDH,  the  contribution  of  horizontal  advection  (the  pro¬ 
cess  of  moving  a  particular  property  of  the  material  from  one  zone 
to  another  by  the  horizontal-velocity  components)  to  the  rate  of 
change  is  calculated.  This  process  depends  only  on  the  horizontal 
components  of  velocity  as  indicated  by  the  following  advection 
equation : 


D tt <p  _  3 tt u <t>  _  3irv<t> 
3t  3x  3y 


(A. 3) 


The  terms  in  Eq.  (A. 3)  are  evaluated  explicitly  in  conserva¬ 
tive  form.  The  difference  equation  is  consistent  with  the  leapfrog 
time  differencing  and  staggered  spatial  grid  of  Figure  (A.2).  In 
finite  difference  notation,  the  advection  terms  are  written: 

■  -  (*x  ^)x  -  (oy  ^)y 

^  -  -  (=y  ~*jx  -  (vx  ^y)y  . 

where  U  =  (u  ,  v)  and  $  =  (0  ,  C).  The  particular  form  of  the 
differencing  is  dictated  by  the  desire  for  a  relatively  simple  con¬ 
servative  form  that  is  accurate  to  second  order  and  is  computationally 


stable.  The  advection  terras  are  evaluated  at  the  midlevel  of 
the  three  time  levels  of  the  leapfrog  scheme,  thereby  assuring 
second-order  time  accuracy. 


Stability  analysis  indicates  that  Eq .  (A. 4)  is  stable  with 
no  damping  when  the  following  inequality  is  observed: 
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Typically  ,  this  stability  condition  is  much  less  restrictive  than 
the  gravity-wave  condition. 

Boundary  conditions  at  the  lateral  mesh  surfaces  are  also 
incorporated  into  the  TENDH  routine.  These  conditions  account  for 
the  ambient  atmospheric  properties  that  enter  the  calculational 
region  from  upwind  and  the  removal  of  the  disturbed  quantities  from 
the  mesh  at  the  downwind  boundary.  These  boundary  conditions  are 
implemented  by  prescribing  ambient  values  on  the  surfaces  that  have 
a  wind  component  directed  into  the  mesh  and  by  extrapolating  values 
from  the  mesh  interior  to  the  boundary  on  those  surfaces  that  have 
wind  components  directed  out  of  the  computational  region. 

PRESUR 

In  the  PRESUR  subroutine,  pressure  gradient  effects  are 
taken  into  account,  and  rates  of  change  of  the  horizontal-velocity 
components  corresponding  to  them  are  formed.  In  the  leapfrog 
equations,  centered  as  indicated  in  Figure  A. 2,  the  pressure  gradients 
can  be  calculated  to  second-order  accuracy  in  both  space  and  time 
without  greatly  complicating  the  equations.  The  calculation  is 
performed  explicitly  (i.e.,  the  gradient  is  evaluated  from  known 
quantities  evaluated  at  the  current  cycle).  Corresponding  to  this 
explicit  formulation,  an  explicit  stability  condition  exists,  which 
severely  limits  the  time  interval.  This  gravity-wave-propagation 
condition  is,  in  fact,  the  most  restrictive  of  all  the  time-interval 


A-8 


constraints  and  typically  limits  the  time  interval  to  a  few  seconds. 
The  stability  condition  is  given  by 


At<  min  (Ax  ,  Ay) 
2/iH 


(A. 6) 


where  H  is  the  effective  thickness  of  the  computational  region, 
and  g  is  the  acceleration  of  gravity. 

The  terms  of  interest  in  the  PRESUR  subroutine  give  rise  to 
a  rate  of  change  of  the  horizontal-velocity  components.  The  differ¬ 
ential  equations  are 
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where  4>  -  gz  =  ^  f  +  gz_  is  the  geopotential,  and  z_  is  the 
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surface  altitude. 

The  above  terms  are  approximated  by  the  difference  equations 
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The  terms  on  the  right-hand  side  of  Eq .  (A. 8)  are  evaluated  from 
quantities  available  at  the  current  time  and  involve  centered  spa¬ 
tial  differences,  which  result  in  second-order  accuracy  in  both 
space  and  time. 

The  calculation  in  PRESUR  is  carried  out  in  two  steps:  first, 
the  geopotential  height  is  formed  by  integration  of  the  hydrostatic 
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equation;  second,  the  contribution  of  pressure  terms  to  the  momen¬ 
tum's  rate  of  change  is  formed.  A  correction  term,  formed  on  the 
first  time  cycle  is  applied  to  each  term  throughout  the  calculation. 
As  seen  in  Eq .  (A.7)fthe  pressure  gradient  is  the  sum  of  several 
terms.  When  terrain  is  present,  each  term  is  large  and  should  can¬ 
cel  each  other  when  the  atmosphere  is  at  rest.  However,  because  of 
the  finite  difference  approximat ion, these  terms  do  not  cancel 
exactly,  giving  rise  to  a  spurious  acceleration.  The  correction 
term  eliminates  this  error,  thereby  increasing  the  accuracy  of  the 
dynamical  calculation. 

RADFLD 

This  subroutine  is  accessed  to  account  for  the  solar  and 
terrestrial  radiation  through  the  atmosphere.  Because  the  amount 
of  time  spent  in  terrestrial  radiation  is  about  one-third  of  the 
total  computational  time,  this  subroutine  is  called  only  at  a  pre¬ 
set  interval  to  update  the  radiation  fluxes.  At  other  times  an 
extrapolation  procedure  is  used  to  approximate  the  radiation 
fluxes . 

RADFLD  calls  subroutine  SOLAR,  which  determines  various 
sun  related  quantities  such  as:  solar  constant  for  current  earth- 
sun  distance,  current  declination  of  the  sun,  current  hour  angle, 
zenith  angle,  normal  insolation, and  other  miscellaneous  quantities. 
RADFLD  then  computes  the  amount  of  water  vapor  in  a  given  vertical 
column,  taking  into  account  the  amount  above  the  column  that  is 
outside  the  mesh  (calculated  using  QEFFG).  Subroutine  RADTER, 
which  determines  the  contribution  to  atmospheric  heating  from 
longwave  radiation,  is  then  called.  The  prescription  used  is 
that  of  Katayama  (1974)  as  included  in  the  3-level  UCLA  GCM.  The 
equations  of  radiative  transfer  are  solved  subject  to  the  follow¬ 
ing  boundary  conditions:  downward  IR  flux  is  zero  at  the  top  of 
the  mesh,  and  the  upward  flux  at  the  earth's  surface  is  blackbody 
at  the  surface  temperature.  The  solution  utilizes  a  bulk  trans¬ 
mission  function  defined  as  the  product  of  water  vapor  and  carbon 
dioxide  transmission  functions.  Absorption  by  ozone  is  neglected 


and  the  transmission  function  is  assumed  linear.  RADSOL  is  then 
called  from  RADFLX,  which  calculates  the  radiant  flux  within  the 
atmosphere  due  to  solar  insolation.  A  one-dimensional  model  is 
used,  which  takes  into  account:  absorption  of  solar  radiation 
by  water  vapor  and  the  Rayleigh  scattering  by  air  molecules. 

TENDV 

As  indicated  in  Figure  A.l,  several  subroutines  occur  within 
a  do-loop  ranging  over  all  of  the  horizontal  zones.  These  sub¬ 
routines  form  quantities  in  which  the  coupling  is  in  the  vertical 
direction  and  are  easily  and  efficiently  calculated  a  column  at 
a  t ime . 

The  first  of  the  calculations  in  the  TENDV  subroutine 
takes  care  of  the  vertical  advection  and  the  rate  of  change  of  the 
surface  pressure.  The  vertical  advection  terms  contain  the  verti¬ 
cal  velocity  (a)  ,  which  is  the  effective  velocity  of  the  fluid 
with  respect  to  the  a-coordinate  .  This  quantity  is  obtained  by 
employing  the  mass  conservation  equation  subsequent  to  the  calcu¬ 
lation  of  the  pressure-tendency  rate. 

The  first  calculation  is  that  of  the  rate  of  change  of 
pressure.  This  quantity  is  obtained  by  summing  the  pressure  equa¬ 
tion  over  all  of  the  vertical  zones.  The  as-yet-unknown-  a  terms 
drop  out  of  the  equation  in  this  procedure 
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where  the  horizontal  velocity  divergence  is  formed  from  quantities 
known  at  the  current  time.  These  terms  are  centered  in  time  and 
space,  so  that  the  pressure  tendency  is  second-order  accurate. 

The  second  calculation  in  the  TENDV  subroutine  is  that  of 
o  .  This  calculation  employs  the  following  equation  in  a  sequen¬ 
tial  fashion. 
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The  first  term  o  (k=l)  is  zero  from  the  top-boundary  condition. 
The  o  values  obtained  in  this  step  are  used  subsequently  to 
evaluate  the  effects  of  vertical  advection  on  all  of  the  primary 
variables.  This  calculation  can  be  considered  to  be  an  explicit 
one;  since  the  velocities  entering  Eqs.  (A. 9)  and  (A. 10)  are  known 
quantities  at  the  current  cycle. 

TENDV  calculates  the  contribution  of  vertical  advection  to 
u  ,  v  ,  6  ,  and  C  using  the  following  relations  respectively: 
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Finally,  Coriolis  contributions  to  the  horizontal  momentum  equations 
are  calculated  in  a  straightforward  manner: 
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The  components  for  the  geostrophic  velocity  ug  and  v^  are 
evaluated  by  interpolation  from  the  ambient  velocity  table. 
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FRICV 


In  preparation  for  the  calculation  of  the  vertical  dif¬ 
fusion  resulting  from  turbulence,  the  boundary-layer  diffusivi- 
ties  are  calculated  in  the  FRICV  subroutine.  They  will  be  used 
subsequently  in  DIFV  to  add  the  effects  of  turbulent  transport 
to  the  tendency  equations.  The  eddy  dif fusivities  for  momentum 
and  heat  are  calculated  according  to  the  ' 2i  level’  prescription 
of  Mellor  and  Yamada  (1974).  In  this  formulation,  a  system  of 
algebraic  equations  involving  the  local  Richardson  number,  tur¬ 
bulent  kinetic  energy,  and  the  local  shear  is  solved  in  order  to 
evaluate  the  eddy  dif fusivities .  The  turbulent  kinetic  energy 
is  described  by  a  differential  equation,  which  amounts  to  solving 
a  tridiagonal  system  of  equations  in  finite  difference  space. 

This  is  done  in  FRICV  by  calling  subroutine  TRIDG.  Kinetic 
energy  is  the  only  quantity  that  is  completely  updated  to 
future-time  level  in  FRICV. 

RADFLX 

The  subroutine  RADFLD,  which  updates  the  radiation  fluxes, 
is  not  called  every  time  step  for  computational  economy;  instead, 
subroutine  RADFLX  is  used  to  update  the  fluxes  by  performing  a 
linear  extrapolation  on  radiative  fluxes  saved  at  two  earlier 
time  steps.  Except  during  sunrise  and  sunset  hours,  this  procedure 
has  been  found  to  work  well. 

DIFV 

This  important  subroutine  serves  several  functions: 

•  The  vertical  turbulent  diffusion  calculation  is  performed 
using  the  diffusion  coefficients  formed  in  the  preceding 
subroutine  FRICV. 

•  The  heat  balance  equation  at  the  soil  surface  is  solved 
efficiently  and  accurately  to  take  into  account  phenomena 
taking  place  on  and  near  the  ground  surface.  An  option 
also  is  provided  for  water/land  surface. 
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As  previously  dicussed,  the  vertical  turbulent  diffusion 
is  a  vigorous  process,  which  acts  across  the  small  zone  thickness 
corresponding  to  the  vertical  zoning.  An  explicit  treatment  of 
this  diffusion  would  require  that  an  additional  numerical  stability 
inequality  be  observed  in  the  determination  of  the  time  interval. 
This  diffusion-stability  criterion  would  be  very  restrictive  and 
would  cause  the  time  interval  to  be  limited  to  even  smaller  values 
than  determined  by  the  other  criteria  discussed  above. 

In  order  to  avoid  such  a  costly  restriction  of  the  time 
interval,  a  numerical  formulation  of  the  terms  involving  vertical 
diffusion  that  is  not  subject  to  the  stability  criterion  was 
implemented.  Such  a  formulation  is  that  of  partial  implicitization , 
in  which  a  linear  system  of  simultaneous  equations  is  solved. 

While  the  method  required  to  evaluate  these  equations  is  somewhat 
more  complicated  (as  described  below)  than  that  for  the  explicit 
formulation,  the  computational  penalty  is  minimal.  The  solution 
is  formed  cellwise  column  by  column  and  makes  use  of  dif fusivities 
calculated  previously  for  the  column  of  cells  in  question.  This 
solution  sequence  dictates  that  the  vertical  calculations  of  the 
FRICV  and  DIFV  subroutines  be  performed  in  the  innermost  loop  of 
the  flowchart  as  indicated  in  Figure  A.l. 

The  terms  of  the  differential  equation  that  are  evaluated 
in  DIFV  are 
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where  <J>  =  (u  ,  v  ,  h  ,  C)  ,  and  <j>  is  a  symbolic  representation 
of  the  other  terms  in  the  governing  equations,  which  have  been 
evaluated  exp’^citly  in  the  subroutines  discussed  above.  The 
difference  formulation  evaluates  the  diffusion  terms  at  the  advanced 
time  of  the  three  levels  of  the  leapfrog  scheme.  This  formulation 
has  been  shown  to  be  unconditionally  stable  against  numerical  error 
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growth.  Taking  into  account  the  variation  in  vertical  cell  size 
with  height,  the  difference  equation  corresponding  to  Eq.  (A. 13) 
is  : 
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In  Eq.  (A. 14),  the  subscript  (k)  indicates  the  vertical  index,  the 
superscript  (n)  indicates  the  time  index,  and  the  subscripts  (+) 
and  (-)  designate  quantities  to  be  evaluated  at  the  interface  of 
the  zone  having  larger  and  smaller  k-indices,  respectively.  All  of 
the  indices  indicating  horizontal  position,  which  are  not  involved 
with  these  terms,  have  been  suppressed. 


Eq .  (A. 14)  is  a  system  of  linear  equations  for  the  unknown 
quantities  +  1  .  which  is  tridiagonal  (the  matrix  of  the 

K  n+1 

coefficients  of  7i<f"  contains  nonzero  elements  only  on  the  main 

diagonal  and  the  bordering  diagonals  immediately  above  and  below 
the  main  diagonal).  This  system  of  equations  is  closed  with  boundary 
conditions  at  the  top  and  bottom  of  the  computational  region.  At 
the  top  of  the  mesh,  zero  turbulent  flux  is  assumed.  At  the  bottom, 
the  velocity  at  the  surface  is  matched  to  a  wall  layer  solution, 
and  the  temperature  boundary  value  at  the  surface  is  assumed  to  be 
known  . 

The  formulation  of  the  surface-heat  balance  in  DIFV  takes 
into  account:  1)  the  sensible  and  latent  heat  fluxes  to  incorporate 
the  logarithmic  dependence  of  wind  speed  on  height,  2)  the  heat  flux 
into  soil,  3)  the  evaporation  and  condensation  of  moisture  from 
the  soil  surlace,  and  4)  simultaneous  solution  of  the  heat  balance, 
atmospheric  temperature  and  water  vapor  equations,  and  the  soil 
temperature  equation.  The  soil  moisture  is  defined  by  a  wetness 
parameter . 


UPDATE 


Subroutine  UPDATE  updates  all  the  primary  variables  to  the 
advanced  time  level.  In  addition,  a  frequency  filter  to  the  time 
integration  (Asselin  1972)  is  applied  at  each  time  step  in  order 
to  remove  the  high-frequency  noise  from  the  solution.  The  filter 
also  helps  to  damp  the  2At  oscillations  inherent  in  leapfrog  time 
differencing  by  averaging  the  solution  over  three  consecutive  time 
steps.  If  the  superscript  (n)  denotes  the  time  level  and  (^) 
denotes  smoothed  variable,  then  the  filter  is  defined  by  the  follow¬ 
ing  operator: 

$n  =  <J>n  +  ia  ($n_1  -  2<|>n  +  <J>n+1)  ♦  (A.  15) 

Note  that  the  filter  is  applied  only  to  old  quantities  ^<f>n) 
after  the  new  quantities  are  calculated.  The  param¬ 

eter  a  =  0  corresponds  to  no  f iltering, and  a  =  1  corre¬ 
sponds  to  a  heavily  biased  filter.  A  value  of  a  =  0.5  com¬ 
pletely  suppresses  2At  oscillations. 

UVBOUN 

Values  of  velocity  at  external  positions  on  the  computa¬ 
tional  mesh  have  not  been  supplied  in  the  previous  subroutines 
and  are  added  in  the  UVBOUN  subroutine.  These  quantities  are 
obtained  in  accordance  with  whether  the  wind  velocity  is  directed 
into  cr  out  of  the  computational  region  at  the  boundary  location 
in  question.  If  the  wind  is  directed  inward,  a  specified  boundary 
value  of  velocity  is  utilized.  On  the  other  hand,  if  the  wind 
is  directed  out  of  the  mesh,  the  boundary  values  of  the  wind 
components  are  obtained  by  extrapolation  from  the  interior  of 
the  mesh. 
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OUTPUT 


The  OUTPUT  subroutine  controls  the  manipulation  of  several 
forms  of  data  that  flow  from  the  computer  to  peripheral  devices. 

This  output  may  take  the  form  of  printer  edits,  tape  dumps,  printer 
plots,  or  film  plots.  The  user  determines  the  frequency  and  mode  of 
data  output  by  specifying  parameters  in  the  INITAL  subroutine. 

These  parameters  specify  how  many  cycles  of  calculation  are  to 
elapse  between  the  successive  edits  of  a  specified  type. 

The  printer  edits  provide  the  user  with  a  record  of  the 
progress  of  the  calculation  in  terms  of  the  values  of  selected 
variables.  These  variables  are  chosen  to  provide  maximum  informa¬ 
tion  about  the  quantities  of  greatest  interest  during  the  normal 
course  of  the  calculation.  More  complete  data  output  is  obtained 
by  requesting  the  so-called  "debug"  edits.  These  edits  contain  a 
large  amount  of  information  on  intermediate  quantities,  which  are 
normally  not  edited.  They  would  be  too  expensive  to  print,  too 
time  consuming  to  inspect,  and  too  cumbersome  to  archive.  The 
debut  edits  find  their  greatest  use  in  diagnosing  abnormal  behavior 
of  a  calculation  by  providing  detailed  information  on  the  suspected 
quantities . 

The  tape  dumps  are  transfers  to  mass  storage  of  all  of 
the  variables  required  to  start  the  calculation.  This  feature 
is  provided  to  enable  the  user  to  reinitiate  a  calculation,  to 
change  data  part  way  through  a  calculation,  or  to  edit  information 
about  a  calculation  selectively  or  in  greater  detail.  These  dumps 
can  be  stored  compactly  for  a  long  period  of  time  and  thus  consti¬ 
tute  a  convenient  and  permanent  record  of  completed  calculations. 

The  printer  plots  provide  the  user  with  various  graphical 
representations  of  the  data.  This  form  of  output  is  invaluable 
in  evaluating  the  results  of  the  calculation.  The  printer  plots 
are  made  during  the  calculation  and  provide  greater  user  overview 
of  the  calculation's  progress  and  insight  into  the  behavior  and 
accuracy  of  the  results. 


The  most  important  and  impressive  form  of  plot  is  on  micro¬ 
fiche.  Contour  plots  of  velocity  magnitude,  temperature,  moisture, 
turbulent  kinetic  energy,  vector  plot  of  the  horizontal  velocity 
vectors  at  different  a  levels,  contour  and  prospective  plots  of 
the  terrain,  etc.,  can  be  generated  at  any  frequency  through  calls 
to  other  plot  routines  from  OUTPUT.  The  plot  routines  are,  however, 
machine  and  computer  center  specific  and  would  require  modifica¬ 
tions  to  utilize  this  capability. 

A. 3  USAGE  AND  INPUT 

As  might  be  imagined,  the  SIGMET  code  that  implements  all 
of  the  methods  and  models  described  above  is  quite  sophisticated 
and  requires  knowledge  and  experience  with  such  techniques  to  oper¬ 
ate  in  a  useful  manner.  It  is  being  continually  updated  in  terms 
of  its  physical  and  numerical  methods  to  enhance  its  usefulness 
and  efficiency.  It,  however,  remains  quite  expensive  to  run, 
presently  requiring  on  the  order  of  10' s  of  minutes  on  a  CDC  7600 
computer  for  a  typical  3-D  mesoscale  simulation  of  approximately 
1-hour  real-time  duration.  It  also  makes  extensive  use  of  com¬ 
puter  storage;  in  its  present  configuration,  it  requires  151. 5gK 
words  to  load  and  130. 5„K  words  of  small  core  memory  to  execute  as 
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well  as  404. OgK  words  of  large  core  memory  (LCM).  In  this  configura¬ 
tion,  it  can  perform  a  calculation  with  a  finite  difference  grid 
design  of  (25  x  21  x  15)  zones  in  the  x,y,  and  z  directions  respec¬ 
tively.  The  sample  calculation  in  the  next  section  provides  coding 
changes  of  array  sizes  required  for  performing  a  slightly  larger 
calculation.  It  is  recommended  that  similar  changes  be  made  for 
a  smaller  calculation  to  minimize  computational  costs.  Because  of 
the  extensive  storage  and  time  requirements,  it  goes  without  saying 
that  conservative  computing  practice  be  followed  when  exercising 
the  SIGMET  code. 

Most  input  variables  required  for  a  SIGMET  calculation  are 
read  through  two  NAMELIST  records:  OUTPT  and  START.  This  input 
uses  standard  CDC  NAMELIST  format.  The  variables  to  be  included  in 
each  record  are  listed  with  a  functional  description  in  Section  A. 5. 


Variables  in  record  OUTPT  control  execution  and  variables  in  START 
provide  numerous  parameters  for  numerical  and  physical  description 
of  a  calculation.  For  an  initial  run,  both  records  should  be  read, 
and  for  a  restart  run  (ISTART  =  1),  when  previously  generated  dump 
tape  is  used  to  initialize  the  calculation,  only  record  OUTPT  is 
read.  SIGMET  automatically  generates  a  dump  tape  on  logical  Unit  8 
for  any  run. 

A  typical  input  deck  for  an  initial  run  is  provided  in  the 
next  section  and  has  the  following  form: 

NAMELIST  OUTPT  :  (see  A. 5) 

NAMELIST  START  :  (see  A. 5) 

UADAT  :  (80  character  field  contain¬ 

ing  upper  air  data  description) 

TRDAT  :  (80  character  field  contain¬ 

ing  terrain  data  description) 

( (HM( I , J ) ,  1=1,  IE),  J=1 ,  JE)  :  terrain  data  in  format 

8F  10.0) 

As  suggested  in  the  above  input  list,  terrain  data  must  be  provided 
for  each  finite  difference  cell  starting  at  the  southwest  corner 
(I,J)  =  (1,1)  for  each  row  in  turn  proceeding  in  an  easterly  (1=1, 

IE)  direction.  These  data  are  read  according  to  the  following  coding 
in  subroutine  INITAL: 

DO  180  J  =  1,  JE 

READ  (5,615)  (HM(I,J),  1=1,  IE) 

180  CONTINUE 

615  FORMAT  (8F10.0) 

As  mentioned  above,  a  restart  run  requires  only  the  input  contained 
in  NAMELIST  OUTPT  to  control  execution  as  well  as  a  dump  tape  of  a 
previous  run  mounted  on  logical  Unit  7. 

A  concrete  example  of  the  input  for  an  initial  run  is  now 
provided  in  the  next  section. 


A.  4  SAMPLE  CALCULATION 

The  setup  and  example  output  for  a  sample  SIGMET  calcula¬ 
tion  are  presented  in  this  section  for  demonstration  purposes.  The 
test  problem  is  the  flow  over  a  three-dimensional  Gaussian  Hill 
with  a  prevailing  uniform  westerly  wind  in  a  neutrally  stratified 
(lapse  rate  =  -0.01  °K/m  =  -  GAM)  atmosphere. 

Table  A.l  below  lists  the  coding  changes  that  are  required 
to  run  this  problem  with  the  present  code  configuration.  These 
changes  include  alteration  in  the  array  dimensions  for  various 
parameters  and  a  code  update  to  calculate  the  hill  terrain  directly, 
rather  than  read  the  terrain  data  from  cards.  All  other  input 
parameters  are  defined  from  NAMELIST  as  described  below.  It  is 
noted  that  the  radiation  arrays  (RADTN,  RADTO,  RADSN,  RADSO)  are 
dimension  to  1,  since  radiation  is  neglected  in  this  test  case 
( IRAD  =  0).  This  results  in  a  substantial  reduction  in  computer 


Table  A.l  Code  Updates  for  Gaussian  Hill  Sample  Case  • 


* IDENT  GAUSS 

•DELETE  COMMON. 2. COMMON.  12 

COMMON/ ARRAY 1/CA! 25.25. 15) .VAI23.25. 15) ,CA(24.24. 15) ,TA<24,24. 15) 
COMMON/ARRAY2/UB125.25, 15) , VB< 25 . 25 . 15) ,CB!24.24. 15) .TB124.24, IS) 
COMMON/ ARRAY3/C! 25. 25. 15) ,V<25.25. 15) ,C( 24.24. 15) .T124.24. 15) 
COMMON/ ARRAY4/CTEN1  25. 25, 15) ,VTEH!25,25, 15) .CTEN(25.25. 15), 

.  TTEN(25.25, 15) 

COMMON/ ARRAY5/PB! (24,24. 15) ,PniPIX!24.24. 15) .PBIPIY124.24, 15) 
COMMON/ ARRAY6/QA! 24,24 , 16) ,QB< 24,24, 16) . HDIFHl 24,24, 15) 
COMMON/ARRAY7/PSA(  24,24) ,PSB<24,24> .PTENC24.24) ,BM(  24,24) , 

PUAl  25,25) , PUB! 25, 25) .PCCT125.25) , CONVARt 25 . 25 ) 
COMMON/ARRAYB/TH(  24,24) ,TC(  24.24) .STEMP! 24,24) ,CC(  24,24) , 

CV(  24,24) , HRATE(24,24) , IJCOLH! 24,24) 

•DELETE  COMMON. 15. COMMON. 18 

COMMON/RAD  I AT/RADTN( 1,1,1) ,RADT0!  1.1,1) ,1UDSN( 1 . 1 . 1) . RADSO ( 1,1,1) 
COMMON/SCRACB/ FACTOR! 24, 16) ,FACT0L!24, 16) ,SDOTR!24, 16) , 

8DOTL!  24, 16) 

•DELETE  COREQV. 4.C0MEQV. 6 

DIMENSION  COMI ! 36036) ,C0M2!36830) , COM3! 36636) .COM4! 37588) , 

.  COM3! 25928) ,C0H6 127072) .C0M714804) .COMB!  4832) . 

.  COM9!  256) . COM 18! 4) , COMI 1!  1536) .C0H12I96) , 

•DELETE  INITAL. 172, IN1TAL. 176 
« ICC- 3868.0 
DO  188  J-l.JE 
DO  188  1*1. IE 
XI- l I-12.5)*DX 
Y1»!J-12.5)*DY 

RH(  I.J)« 1806.0*EXP(-!XI»«24YI«*2)/8ICC**2> 

188  CONTI NOE 
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core  storage  requirements.  Further  relief  could  be  obtained  by 
deleting  the  plotting  subroutines  (CONTUR,  PLOTER,  AND  VELVER)  and 
adding  the  dummy  subroutines  MODESG  and  EXITG ,  since  microfiche 
plots  are  not  to  be  generated  (IFILM  =  0).  This  step  should  be 
implemented  if  SIGMET  is  to  be  run  at  a  computer  facility  where  the 
identical  plotting  package  is  not  available. 

The  following  are  the  NAMELIST  input  parameters  and  descrip¬ 
tion  caption  input  for  this  sample  calculation  : 


$OUTPT 

NCYCLS=5 ,  IEDIT=5 ,  IOUT=0,  IPLOT=0,  ISTART=0, 

KBEG=15 ,  KBGPL=15 ,  IFILM=0, 

$END 

$START 

GMTMHR=7 . 25 ,  KMAX=15,  DELTAT=1 . 6 ,  DLAT=0.0, 

ANGLE=0 . 0 ,  DLONG=106 . 4 ,  PTOP=650.0,  PBOT=1000., 

W=0.0,  IRAD=0 ,  XORG=0 . 0 ,  XCLAY=0.0, 

XSAND=1 . 0 ,  RPDIF=0 . 608 ,  CPDIF=2.0,  TRANS=1.0, 

VEG=0 . 0 ,  ZNAUT=0 . 1 ,  QEFFG=0 . 5 ,  GDEPTH=0.1, 

TDEEP=15 . 0 ,  TWATER=15 . 0 ,  ZBL=200.0,  SUREM=1.0, 

IE=24  ,  JE=24  ,  DX=750 . 0 ,  DY=750.0, 

TG=576*15 . 0 ,  GW=576*0 . 2 ,  HRATE=576*0 . 0 ,  I JCOLM=576*0 , 

GAM=0 . 01 ,  GAM1=0 . 005 ,  PINV=500.0,  TNPRS=15.0 

NPRES=21 ,  DAYBEG=1 . 0 ,  ZBOT=0.0,  ALPH=0.0, 

PPRES=600 . ,620. ,640. ,660. ,680. ,700. ,720. ,740. ,760. ,780. ,800. , 

820. ,840. ,860. .880. ,900. ,920. ,940. ,960. ,980. ,1000. , 

UPRES=21*10. 0 ,  VP RES=2 1*0.0,  CPRES=21*0 . 0 ,  THETA=21*270 . , 

RELHM=21*0 . 0 ,  CX3=0 . 013882 ,  CX4=1.2, 

$END 

SIMULATED  INITIAL  CONDITIONS 
3-D  GAUSSIAN  HILL 

If  this  were  a  usual  SIGMET  run  in  which  terrain  data  were  to  be 
input  via  cards,  the  above  input  would  be  followed  by  terrain  data 
cards  in  an  8F  10.0  format. 

The  above  completes  the  coding  changes  and  data  input  for 
the  sample  case.  Table  A. 2  is  a  reproduction  of  the  raw  computer 
output  for  this  calculation.  The  numbers  and  tables  provided  there 
are  relatively  self  explanatory.  For  illustration  and  check-out 
purposes,  only  the  first  two  pages  of  output  for  the  initial  condi¬ 
tions  and  for  the  results  at  the  end  of  the  5-cycle  run  are  provided. 


Table  A. 2  Sample  Output  (Initial  Condition), 

(cont inued ) 
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Table  A. 2  Sample  Output  (5-cycle  run), 

(continued) 
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GLOSSARY  OF  INPUT  VARIABLES 


A. 5 

The  set  of  input  variables  required  to  define  a  SIGMET 
simulation  is  presented  in  the  following  table.  Variables  are 
input  in  standard  CDC  namelist  format  in  either  NAMELIST  /OUTPT/ 
or  NAMELIST  /START/  as  indicated. 

NAMELIST  /OUTPT/ 

Variable  Description 

ISTART:  Starting  parameter 

=  0  for  initial  run 

■»  0  for  restart  runs.  For  restart  runs  save  tape 
should  be  mounted  on  logical  Unit  7 


IEDIT: 

Frequency  of 

edit  cycle 

I PLOT: 

Frequency  of 

printer  plot 

I  OUT : 

Frequency  of 

debug  edit  cycle 

I FILM: 

Frequency  of 

film  plot 

NCYCLS : 

Total  number 

of  time  cycles 

KB  EG : 

Starting  edit 
(if  IEDIT  is 
by  level  for 

level 

nonzero  printer  edit  is  provided  level 
all  variables  from  K  =  KBEG  to  KMAX) 

KBGPL : 

If  IFILM  is  nonzero  plots  on  horizontal  slices 
K=KBGPL  to  KMAX  are  to  be  prepared 

NAMELIST  /START/ 

Variable  Description 

KMAX:  No.  of  vertical  levels 

IE:  No.  of  computational  zones  in  x-direction 

JE:  No.  of  computational  zones  in  y-direction 

DX:  Grid  size  in  x-direction  in  meters 

DY:  Grid  size  in  y-direction  in  meters 

DF.LTAT:  Time  step  size  in  seconds 
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Variable 


■* 


GAM 

TNPRS 

GAM1 

PINV 


Description 

By  assigning  a  non-negative  value  to  GAM  ,  the  measured 
profile  of  temperature  (TPRES)  can  be  replaced  by  one 
with  TNPRS  (in  °C)  at  ground  decreasing  with  altitude 
at  lapse  rate  GAM  (in  °C  per  meter).  An  inversion 
layer  can  be  simulated  by  having  a  smaller  lapse  rate 
GAM1  (in  °C/meter)  starting  at  a  level  where  pressure 
is  PINV  (in  millibar). 


ANGLE : 

ZBOT: 

TG: 


GW: 


HRATE : 

TWATER : 

TDEEP : 

XORG  I 

XCLAY 

XSANDI 

VEG: 

ZNAUT: 

DLAT : 
DLONG : 

Q. 

CXI,  CX2 : 
CX3 ,  CX4 : 


x-axis  rotation  measured  counter  clockwise  from  East 
in  degrees 

Ground  elevation  (AGL)  of  measurement  location  if 
available 

A  two-dimensional  array  at  least  IE  x  JE  long  repre¬ 
senting  ground  temperature  is  °C.  If  measured  values 
for  TG  are  not  available,  they  can  be  extrapolated 
from  the  measurement  of  temperature  in  the  atmosphere. 
This  is  done  by  setting  the  value  of  GDEPTH  negative. 

A  two  dimensional  array  at  least  IE  x  JE  long  repre¬ 
senting  ground  moisture  expressed  as  fraction  of  the 
saturated  humidity  of  air  near  the  ground 

A  two-dimensional  array  at  least  IE  x  JE  long  repre¬ 
senting  external  ground  heat  source,  if  any 

Temperature  in  °K  of  water  body,  if  any,  present  in 
the  region  of  simulation 

Temperature  in  °K  of  underground  soil  assumed  constant 
everywhere 

Composition  of  soil  in  the  region  of  simulation,  in 
fraction  by  mass  of  organic  materials,  clay  and  sand, 
respectively 

Vegetated  fraction  of  the  total  surface  area 

Roughness  length  in  meters  typical  to  the  terrain  over 
the  region 

Latitude  in  degrees  of  region  of  simulation 

Longitude  in  degrees  of  region  of  simulation 

One-dimensional  array  of  at  least  KMAX  +1  long  for 
o  values 

Two  parameters  specifying  the  log-linear  distribution 
of  o  levels 

Two  parameters  specifying  o  values  according  to  a 
geometric  progression  series 

Note:  If  one  of  two  parameters  is  zero,  o  values 

are  not  computed  and  preset  values  (read 
through  Q)  are  used. 
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Variable 


GMTMHR 

DAYBEG : 
PTOP: 

PBOT : 

UPRES 
VP  RES 
TP  RES 
CP  RES 

PPRES : 

NPRES  : 
THETA : 


RELHM: 


ALPH: 
ZBL : 


I  RAD: 

(W  ) 

)  TRANS l 

) SUREM ( 
( QEFFG  ) 

GDEPTH : 


Descript  ion 

Gmt  hour  (in  fraction  notation)  corresponding  to 
initial  simulation  time 

Julian  day  of  initial  simulation 

Pressure  in  millibar  corresponding  to  the  top  level 
of  the  model 

Standard  sea  level  pressure  in  millibar 
Ambient  (usually  obtained  from  soundings)  values  of 
horizontal  wind  in  meter/sec  or  in  mi/h,  temperature 
°C,  humidity  in  kg/kg.  These  are  one-dimensional 
arrays  of  length  30  each. 

One-dimensional  array  of  size  30  representing  pressure 
in  millibar  where  UPRES,  VPRES, etc . ,  are  measured 

No.  of  entries  of  UPRES,  VPRES, etc.  (<.  30) 

A  one-dimensional  array  of  at  least  NPRES  length 
representing  angle  in  degrees  of  wind  vector  measured 
clockwise  from  north.  If  the  measurement  is  available 
in  magnitude/direction  format,  UPRES  should  contain 
the  magnitude  in  meters  per  sec  and  THETA  the  direction. 
If  component  measurements  are  to  be  entered,  UPRES 
should  contain  the  East-West  component ,  and  VPRES  should 
contain  the  North  South  component,  both  in  miles  per 
hour,  and  all  THETA  should  be  set  at  zero. 

An  array  of  at  least  NPRES  elements  representing  values 
of  relative  humidity  at  NPRES  levels  of  the  atmosphere. 
Expressed  as  percentage  of  saturated  value  for  air  at 
given  temperature.  If  any  of  RELHM  is  nonzero,  the 
measurement  CPRES  is  replaced  by  new  values  calculated 
from  RELHM  and  TP RES. 

Time  smoothing  parameter  to  damp  out  high  frequency 
temporal  oscillations,  01ALPH11 

Boundary  layer  height  in  meters  above  which  an  inter¬ 
polated  value  from  UPRES,  VPRES  is  used  for  initiali¬ 
zing  the  velocity  field.  Below  ZBL  a  power  law  profile 
is  used. 

Frequency  of  cycles  when  radiation  subroutines  are 
called.  If  zero,  radiation  contribution  is  completely 
ignored . 

Radiation  parameters 


Length  scale  in  meters  required  for  solution  of  ground 
temperature 
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Variable 


Description 


i 


IJCOLM 


A  two-dimensional  array  of  size  IE  x  JE,  each  element 
corresponding  to  a  given  vertical  column.  A  nonzero 
value  will  result  in  printout  of  a  column  (liKiKMAX) 
of  flow  variables  every  IOUT  cycles. 
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APPENDIX  B 
UNMET  ORGANIZATION 

B.l  LINMET  ORGANIZATION 

The  3-D  steady  linear  code  LINMET,  to  calculate  strati¬ 
fied  flow  over  terrain  was  developed  during  the  contract.  This 
code,  in  essence,  solves  the  Scorer’s  equation  for  the  Fourier 
transform  of  the  perturbation  vertical  velocity  and  then  computes 
the  perturbation  to  all  three  components  of  the  velocity.  As 
indicated  in  the  flow  diagram  (Table  B.l),  the  various  specific 
tasks  are  done  through  calls  to  different  subroutines  from  the 
main  program.  The  code  is  written  in  FORTRAN  IV  ensuring  porta¬ 
bility  to  any  computer.  An  exception  is  the  DPLT  subroutine  that 
makes  use  of  machine  dependent  plotting  routines  available  at 
Lawrence  Berkeley  computing  facility.  This  is  the  only  subroutine 
that  has  to  be  deleted  or  replaced  for  adaptation  to  other  machines. 

B. 2  SUBROUTINE  DESCRIPTION 

In  the  following, we  give  a  short  description  of  each  of  the 
subroutines  in  LINMET. 

INPUT  SUBROUTINE 

The  INPUT  subroutine  processes  data  provided  through  NAMELIST 
data  cards,  so  that  all  information  required  for  subsequent  normal 
execution  of  a  problem  simulation  are  available.  These  data  cards 
provide  values  for  the  following  quantities: 

•  spatial  step  sizes  both  in  horizontal  and  vertical 
direction 

•  powers  of  two,  which  define  the  number  of  grids  in 
the  horizontal  directions 

•  the  number  of  grids  in  the  vertical  direction 

•  height  of  the  computing  mesh 

'  •  editing  and  plotting  parameters 

•  option  switch  for  self-prescribed  standard  data 

The  terrain  height  data  are  normally  read  from  a  disk  file.  INPUT 
selects  a  limited  area  from  the  entire  terrain  data  and  then  per¬ 
forms  a  boundary  smoothing  to  remove  boundary  related  errors 
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TABLE  B.l  LINMET  Logical  Flow  Diagram. 
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arising  from  periodic  boundary  condition.  It  then  computes  the 
Fourier  transform  of  the  smoothed  terrain  data  for  later  use. 

It  also  computes  the  constant  vertical  mesh  with  the  average 
terrain  height  as  the  lowest  level  and  the  top  of  the  computing 
mesh  as  the  highest  level.  Through  a  switch,  INPUT  is  capable 
of  providing  sel f -prescribed  meteorological  and  terrain  data 
that  can  be  used  for  checking  out  various  aspects  of  the  code. 

INTERR  SUBROUTINE 

The  meteorological  data  read  in  subroutine  INPUT  are  .  ot 
necessarily  at  the  levels  of  the  vertical  mesh.  INTERR  inter¬ 
polates  these  data  to  the  levels  of  the  vertical  grid  and  stores 
them  in  arrays.  It  also  computes  the  derived  quantities  from 
these  data  arrays  and  stores  them  for  later  use. 

SCORER  SUBROUTINE 

This  subroutine  is  called  from  the  main  program  in  a  loop 
covering  all  the  wave  numbers.  SCORER  first  computes  the  multi¬ 
pliers  in  the  Fourier  plane  ,  which  are  related  to  various  finite 
difference  operations  in  the  physical  plane.  It  then  computes  the 
Scorer  parameter  and  the  wave  vector  and  stores  them  in  arrays. 

In  doing  so  it  checks  for  critical  levels,  if  any,  and  thereby 
sets  up  the  domain  of  integration  (from  ground  to  critical  level, 
if  any,  or  top  of  the  computing  region).  It  also  sets  up  a  switch 
to  avoid  integration  for  certain  modes  for  which  the  solution  iv 
trivial. 

C.'ALW  SUBROUTINE 

This  subroutine  controls  the  other  subroutines  that  inte¬ 
grate  the  Scorer's  equation  to  solve  for  the  perturbation  vertical 
velocity.  It  first  distinguishes  between  two  categories  based  on 
the  behavior  of  the  Scorer's  equation  at  the  top  boundary.  It 
then  sets  up  the  initial  values  of  the  functions  corresponding  to 
the  top  boundary  before  calling  the  integration  routines.  The 


integrated  function  and  its  derivatives  are  stored  for  each  vertical 
level.  The  latter  is  required  to  compute  the  horizontal  perturba¬ 
tion  velocity  from  the  vertical  perturbation  velocity.  Once  the 
integration  is  complete,  the  bottom  boundary  condition  is  utilized 
to  give  the  actual  vertical  perturbation  velocity  in  the  Fourier 
plane.  Before  returning  to  the  main  routine,  CALW  calculates  the 
other  perturbation  components  as  well. 

INT  SUBROUTINE 

Subroutine  INT,  which  is  called  from  CALW,  performs  a  step- 
by-step  solution  of  a  system  of  first-order  differential  equations 
(initial  values  problem)  by  sixth-order  Adams-Moulton  Predictor 
Corrector  method,  with  starting  procedure  based  on  Rosser's  formula. 
It  employs  double  precision  arithmetic  and  uses  variable  step  size 
to  satisfy  a  preset  error  criterion,  but  steps  are  forced  to  land 
on  multiples  of  a  specified  print  interval.  Before  INj  is  called, 
another  subroutine  (.INTO)  must  be  called,  which  initializes  various 
constants  used  by  INT. 

DERIVA,  DERIVB  SUBROUTINES 

These  subroutines  are  called  from  INT  to  provide  derivatives 
of  the  function  at  given  vertical  height.  DERIVA  subroutine  supplies 
the  derivatives  for  Case  A  integration,  and  DERIVB  does  the  same  for 
Case  B.  These  subroutines  are  declared  in  EXTERNAL  statements  in 
CALW,  because  their  names  appear  as  dummy  parameters  in  calls  to  INT 
routine. 

CAPPZ  SUBROUTINE 

This  subroutine  function  provides  the  value  of  the  modified 
wave  vector  at  any  altitude  by  interpolation  from  the  wave  vector 
table  that  was  set  up  by  SCORER.  This  is  frequently  called  from 
DERIVA/DERIVB  subroutines  that  must  supply  derivatives  to  INT 
at  a  great  many  subdivisions  within  a  vertical  grid  to  satisfy  the 
error  criterion. 
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SCORZ  SUBROUTINE 


This  subroutine  function  is  identical  in  structure  to  CAPPZ . 
It  supplies  interpolated  values  of  the  Scorer's  parameter  at  any 
altitude  to  DER IVA/DERIVB . 

REFFT  SUBROUTINE 

Once  the  perturbation  components  are  calculated  for  all  wave 
numbers  and  for  all  vertical  levels,  this  subroutine  is  called  from 
the  main  program.  REFFT  calls  F0UR2  at  each  vertical  level  to  per¬ 
form  the  inverse  Fourier  transform.  The  results  are  stored  in 
three-dimensional  arrays. 

F0UR2  SUBROUTINE 

This  is  a  fast  and  efficient  Fourier  transform  routine 
employing  Cooley-Tukey  algorithm  and  capable  of  handling  multi¬ 
dimensional  data.  The  dimension  of  the  data  is  required  to  be  an 
integral  power  of  two.  Switches  are  provided  so  that  the  same 
routine  can  handle  real,  complex,  and  half-complex  data  and  do  both 
Fourier  synthesis  and  analysis. 

OUTPUT  SUBROUTINE 

The  OUTPUT  subroutine  controls  the  manipulation  of  two  forms 
of  data  that  flow  from  the  computer  to  peripheral  devices.  These 
outputs  are  in  the  form  of  printer  edits  or  film  plots.  Both  are 
activiated  through  switches  read  from  data  cards.  Printer  edit 
provides  the  three  perturbation  velocity  components  along  each  hori¬ 
zontal  slice  starting  from  the  top.  Film  plots  include  fcr  each 
horizontal  slice  a  vector  plot  and  two  contour  plots:  one  for  the 
magnitude  of  the  horizontal  velocity  vector  and  the  other  for  ver¬ 
tical  velocity  component.  On  the  background  of  each  frame,  the  con¬ 
tour  of  the  terrain  is  provided. 


DPLT  SUBROUTINE 


This  subroutine  prepares  the  graphic  output  on  film  by 
making  use  of  various  plotting  subroutines  available  through  IDDS 
graphic  package  at  Lawrence  Berkeley  Laboratory.  In  the  event  of 
adaptation  of  this  code  to  a  different  computer,  the  contents  of 
this  subroutine  should  be  thoroughly  reviewed,  and  changes  should 
be  made  to  utilize  the  available  plotting  routines.  If  plots  are 
not  required,  the  whole  subroutine  can  be  replaced  by  a  dummy 
subroutine  by  the  same  name. 

B . 3  A  SAMPLE  LINMET  RUN 

The  deck  setup  and  the  output  from  a  typical  LINMET  run  are 
included  here  for  demonstration.  The  test  problem  over  a  grid  of 
32  x  4  x  10  is  initialized  by  self-prescribed  standard  data  pro¬ 
vided  in  INPUT  (by  setting  IDUMDAT=1  in  NAMELIST  data).  This  self- 
prescribed  data  consist  of  terrain  heights  corresponding  to  a 
two-dimensional  Gaussian  mountain  (no  variation  on  y-direction) ,  a 
uniform  westerly  wind  of  10  meters/second,  and  a  linear  temperature 
profile  corresponding  to  a  constant  lapse  rate  with  T  =  288°K  at 
sea  level . 

Table  B.2  lists  the  change  cards  required  to  modify  the 
existing  version  of  the  code  for  this  test  run.  Notice  that  the 
references  to  all  plotting  routines  in  DPLT  are  deleted,  since  plots 
were  not  requested  for  this  run.  Table  B.3  lists  the  input  data 
cards.  Notice  that  with  the  option  of  IDUMDAT=1 ,  the  meteorological 
data  are  not  required.  Following  Table  B.3  is  the  reproduction  of 
the  raw  computer  output  for  this  run. 
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Table  B.2  List  of  Change  Cards  Required  for  the  Demonstration  Run 


HUNT  TEST 

*  I  DENT  TEST  _ _  _  _  _ 

.  •IiKI.KTK  COMMON :  15. COMMON  .  1 H 

COMMON /TKRA 1  N/T<  )PO  (4  1  ,  5 1  )  ,il  (32  ,  4 ) 

COMMON /F 1  K1JJK /U  (34 ,4 .20)  ,V  (34 ,4 ,20)  ,  W(34  .4  ,20) 

COMMON /KKTAl.L/UTRANZ  (2,17,4)  .  VTRANZ  (2,17,4)  ,  WTRAN7,(2 ,17,4), 
ZTRANZ  (2,17,4)  ,N  (2) 


•OKI. KIT:  COMMON.  20 

DIMKN.3  ION  REFFTU  (32,4)  , REFFTV  (32 ,4)  . REFFTW(  32,4  )  ,1111(32 ,4) 
^DELETE  DPI.T.4  ,DrLT.218  _ _  _ _ _ 


Table  B.3  NAMELIST  Input  Data  for  the  Demonstration  Run. 


$START 

DX 

= 

. 1 E+04  , 

DY 

= 

. 1E+04 , 

MPOW 

= 

5, 

NPOW 

- 

2, 

NZ 

= 

20, 

ZTOP 

= 

. 5E+04  , 

NSTN 

= 

1  , 

MBEGPR 

= 

18  , 

MBEGPL 

= 

18, 

IEDIT 

= 

1 . 

I  FILM 

= 

0. 

IDUMDAT 

= 

1 . 

SEND 
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Copy  available 

Penr.ii  f-ljy  p.,7 


to  DTIC  does  rwj 
•s!c  Joproductiei. 


Table  B.4  Computer  Output  for  the  Demonstration  Run. 
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Table  B.4  (cont’d)  Computer  Output  for  the  Demonstration  Run. 
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VARMET  ORGANIZATION 


C.l  VARMET  ORGANIZATION 

VARMET  is  written  in  modular  form  to  enhance  the  reading 
of  the  code  and  to  facilitate  structural  changes  or  additions. 

The  code  is  divided  into  three  phases  as  shown  in  Figure  C.l. 

The  initialization  phase  defines  the  problem  parameters,  generates 
the  computational  mesh,  and  specifies  the  estimated  flow  field. 

The  adjustment  phase  computes  a  nondivergent  flow  field  based  on 
the  estimated  field  determined  in  the  initialization  phase,  and  the 
the  output  phase  edits  the  flow  field  variables  and  saves  the  final 
solution  on  tape.  Figure  C.2  is  a  flowchart  showing  the  logical 
flow  of  the  Fortran  code  in  addition  to  summarizing  the  subroutine 
f unct ions . 

C.2  SUBROUTINE  DESCRIPTIONS 

CARDS  SUBROUTINE 

CARDS  interprets  free-form  input  card  data,  loading  blank 
common  according  to  the  position  of  a  variable  within  the  common 
block.  CARDS  performs  a  function  similar  to  a  namelist  statement, 
but  unlike  namelist,  it  is  machine  independent.  Section  C.3  describes 
the  CARDS  input  form  accepted  by  subroutine  CARDS. 

CNWAD  SUBROUTINE 

CNWAD  adjusts  the  estimated  flow  field  in  conformal  space. 

The  estimated  wind  components  are  first  transformed  from  real  to 
conformal  space.  The  velocity  potential  is  then  determined  based 
on  the  transformed  wind  components  through  successive  line  over¬ 
relaxation  algorithm.  Once  the  potential  field  is  known  the  non¬ 
divergent  velocity  field  is  computed.  Finally,  an  inverse  trans¬ 
formation  is  applied  to  determine  the  final  values  of  the  velocity 
components  in  real  space. 
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START 


Sec  Default 
Problem  and 
Program  Con¬ 
trol  Parameters 


▼  INPUT 

Read  Problem 
and  Program 
Control 
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SETUP 
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Finite 

Difference 

Mesh 


▼  CARDS 

Read  Free 
Form  Data 
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Print  Adjusted 
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DBUG-1 
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Windfield  in 
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Windfield 


FINISH 


Solve  for 
Adjusted 
Windfield 
Via  SOR 
Itera 


Figure  C.2  Subroutine  Flow  Diagram. 
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DEFINE  SUBROUTINE 


DEFINE  zeroes  common  area  and  sets  default  values  of  pro¬ 
gram  and  problem  control  parameters.  A  list  of  default  values  is 
given  in  C.4. 

EDIT  SUBROUTINE 

EDIT  generates  line  printer  output  of  the  flow  field 
variables.  A  control  flag  may  be  set  to  bypass  this  capability 
if  the  data  are  processed  by  a  separate  program  developed  to  pro¬ 
duce  vector  and  contour  plots  of  the  calculated  wind  field  or  if 
output  produced  by  subroutine  PEDIT  is  preferred.  EDIT  may  be 
called  from  two  locations  in  the  program  flow  (see  Figure  C.2). 

EDIT  prints  for  each  cell  (i,j,k)  the  velocity  components  (u,v,w), 
the  local  divergence,  and  velocity  potential  in  blocks  correspond¬ 
ing  to  horizontal  planes.  All  or  part  of  the  computational  domain 
may  be  printed  by  appropriately  setting  IMN,  IMX,  JMN,  JMX,  KMN, 
and  KMX. 

ENDLCH  SUBROUTINE 

ENDLCH  adjusts  the  estimated  flow  field  in  real  space 
through  an  SOR  iteration  to  obtain  the  divergence-free  wind  field. 
The  effects  of  topography  and,  to  some  extent,  atmospheric  stability 
are  accounted  for  in  the  iteration.  Boundary  conditions  are  imposed 
only  along  the  terrain  boundary  where  the  normal  velocity  is  taken 
to  be  zero.  Velocities  at  the  lateral  and  upper  boundary  surfaces 
are  determined  by  the  iteration  procedure. 

INPUT  SUBROUTINE 

Terrain,  surface  winds,  mesh  parameters,  and  program  control 
data  are  entered  through  subroutine  INPUT.  INPUT  utilizes  sub¬ 
routine  CARDS  to  initialize  the  mesh  parameter,  program  control, 
and  surface  winds  data.  Terrain  data  may  be  entered  through  an 
input  tape,  data  cards,  or  data  statement. 
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IinTERP  subroutine 


INTERP  generates  the  horizontal  distribution  of  the  estimated 
wind  field.  The  observed  data  at  a  given  level  above  the  terrain 
are  interpolated  using  al/r2  weighting  factor  to  obtain  the  set  of 
velocity  components  corresponding  to  the  appropriate  cell  faces  for 
that  level. 

MAP  SUBROUTINE 

MAP  generates  line  printer  contour  plots  of  the  terrain, 
velocity  components,  wind  speed,  divergence,  and  velocity  potential. 
For  three-dimensional  domains,  contour  plots  are  defined  for  the 
distribution  of  a  variable  on  a  horizontal  plane.  When  the  two- 
dimensional  option  is  chosen,  the  contour  plots  represent  the  dis¬ 
tribution  in  the  2-D  plane  of  the  computation. 

OBSWND  SUBROUTINE 

OBSWND  initializes  the  estimated  flow  field.  Surface  data 
from  one  or  more  observation  stations  and/or  upper-air  data  from  one 
observation  station  are  interpolated  to  define  the  discrete  field 
of  estimated  velocity  components  throughout  the  mesh. 

PEP IT  SUBROUTINE 

PEDIT  provides  a  summary  edit  of  mesh  parameters  and  calls 
subroutine  MAP  to  provide  contour  line  printer  plots  of  the  flow 
variables.  PEDIT  may  be  used  in  place  of  subroutine  EDIT  to  reduce 
the  volume  of  output  or  if  detailed  numbers  are  not  required. 

SETUP  SUBROUTINE 

Subroutine  SETUP  initializes  the  finite  difference  mesh 
according  to  mesh  spacing  and  domain  size  information  provided  by 
INPUT.  SETUP  assumes  constant  zoning  in  both  the  x  and  y  directions. 
For  real  space  calculations,  zoning  in  the  vertical  direction  is 
also  assumed  constant.  In  the  conformal  mode,  either  constant  or 


variable  zoning  may  be  used.  Variable  zoning  may  be  generated  in 
one  of  two  ways.  A  log-linear  distribution  is  available  that 
assumes  a  minimum  zone  size  at  the  surface  and  expands  the  zoning 
with  increasing  altitude  according  to  a  logarithmic  formula.  In 
lieu  of  this  expansion,  a  geometric  increase  in  zone  size  may  be 
chosen.  Again,  the  model  assumes  a  minimum  zone  size  at  the  terrain 
surface  and  expands  subsequent  zones  in  the  vertical  by  a  constant 
factor,  i  computed  in  SETUP  based  on  DZMIN,  KMAX  ,  and  ZDIST.  To 
avoid  large  truncation  errors,  <  should  be  kept  in  the  range 
1 . 0<  y  <1 . 3 . 

TAYLOR  SUBROUTINE 

TAYLOR  interpolates  a  wind  speed  and  direction  at  a  speci¬ 
fied  number  of  points  within  the  mesh  using  a  second-order  TAYLOR 
expansion  about  the  points  of  interest.  TAYLOR  has  no  effect  on 
the  solution  of  the  adjusted  wind  field  and  is  only  used  for 
editing  purposes. 

TRIDG  SUBROUTINE 

TRIDG  inverts  the  tridiagonal  matrix  composed  of  the  coeffi¬ 
cients  defined  by  the  finite  difference  approximation  to  the  govern¬ 
ing  equation  for  the  velocity  potential  along  a  column  (ij).  TRIDG 
is  a  special  adaptation  of  the  Gaussian  Elimination  Procedure. 

WTAPE  SUBROUTINE 

Subroutine  WTAPE  dumps  the  common  block  onto  tape  after  the 
estimated  and  adjusted  wind  fields  are  calculated.  The  common  block 
contains  all  the  pertinent  zoning  information  and  field  variable 
arrays . 
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DESCRIPTION  OF  INPUT  DATA 


VARMET  requires  several  types  of  input  data  to  initialize 
a  run.  Control  flags  are  used  to  guide  the  simulation  along  one 
of  many  ;aths  present  in  the  code.  Mesh  parameters  are  necessary 
to  describe  the  problem  geometry.  Terrain  data  are  needed  to 
define  the  bottom  boundary  of  the  computational  domain,  and 
observed  wind  data  are  required  to  determine  an  estimated  wind  field 
throughout  the  domain. 

Free-form  input-card  data  comprising  the  control  flags, 
mesh  parameters,  and  surface/upper  air  wind  observations  are  read 
by  subroutine  CARDS.  Hollerith  headings  for  specified  locations 
within  the  mesh  where  wind  speed  and  direction  data  are  to  be 
edited  (see  below)  are  read  from  cards  during  the  initialization 
phase.  Terrain  data  are  read  from  cards  or  tape  or  initialized 
via  a  data  statement  in  subroutine  INPUT.  The  format  structure 
of  terrain  data  may  need  to  be  modified  in  INPUT  to  conform  to  the 
specific  data  tapes  being  processed. 

The  Input  Variable  Glossary  in  Section  C.4  contains  the 
necessary  input  parameters  read  by  subroutine  CARDS.  These  variables 
comprise  the  start  or  initial  simulation  deck.  VARMET  has  the  capa¬ 
bility  to  process  one  or  more  wind  simulation  cases  in  a  single  run. 
If  operating  in  this  mode, additional  blocks  of  free-form  data  are 
included  behind  the  start  deck  and  serve  to  modify  only  those  param¬ 
eters  that  change  from  simulation  to  simulation.  A  detailed  descrip¬ 
tion  of  the  general  card  input  form  and  data  block  arrangement  is 
described  below. 

FORMAT  OF  THE  INPUT  DATA  CARDS 

The  input  data  cards  found  in  the  glossary  are  read  by  sub¬ 
routine  CARDS  according  to  the  following  format: 

1.  An  asterisk  is  entered  in  column  1. 

2.  After  the  asterisk  any  form  of  a  Hollerith  label 
(normally  the  variable  name)  may  be  entered. 

3.  A  dollar  sign  is  entered. 


t 


The  location  in  Blank  Common  is  entered 

followed  by  a  comma. 

The  numerical  value  of  the  variable  is  entered; 

a.  if  the  variable  is  integer,  enter  the  value 
as  integer. 

b.  if  the  variable  is  real,  enter  the  value  in 
floating  point  or  E-format. 

c.  if  the  variable  is  logical,  enter  a  1  for 
true,  0  for  false. 

d.  if  an  array  is  being  input,  either  enter 
each  element  individually  (as  in  a-c  above), 
or  enter  multiple  elements  on  one  card. 


Example : 


f  *WDR  $41,320.0,270.0,340.0 
♦ORIGNX  $35~ 3. 521+05  — 


[  »I2D  $23,0 
(  *DX  $6,2500,0 
•BETA  $3,0.14 


To  end  a  block  of  data,  the  field  following  the  dollar  sign  is 
left  blank.  For  example: 


•END  OF  START  BLOCK  $ 


r 


The  start  deck  should  contain  values  for  each  of  the  vari¬ 
ables  listed  in  the  glossary  that  differ  from  their  preassigned 
value  specified  in  subroutine  DEFINE.  If  more  than  one  simulation 
will  be  processed,  additional  decks  for  each  subsequent  simulation 
must  be  included.  These  additional  decks  should  contain  only  those 
variables  that  need  modifying.  In  most  cases  the  modifications 
will  include  at  a  minimum  the  redefinition  of  the  date  and  time  of 
the  simulation  and  of  the  observed  wind  station  data.  An  example 
of  a  multisimulation  deck  setup  is  given  in  Figure  C.3. 

It  is  often  desirable  to  edit  wind  speed  and  direction 
data  at  specified  locations  within  the  computational  domain  that 
do  not  necessarily  coincide  with  a  mesh  boundary.  This  capability 
is  provided  by  subroutine  TAYLOR.  TAYLOR  computes  and  edits  wind 
information  along  a  vertical  column  for  up  to  25  specific  locations 
within  the  area  of  interest.  Hollerith  headings  (of  up  to  four 
characters)  may  be  provided  via  data  cards  for  these  locations. 
These  cards  must  appear  after  the  first  simulation  deck  and  before 
the  terrain  data  deck  as  shown  in  Figure  C.3.  This  Hollerith 
information  must  be  entered  in  a  20A4  format.  Default  numerical 
headings  will  be  provided  if  these  cards  are  not  present. 

Terrain  data  (assumed  to  be  in  meters)  may  be  initialized 
in  one  of  three  ways  as  discussed  above.  Data  on  tape  are  read 
from  unit  8  (8F10.8  format).  Data  on  cards  are  read  from  a  ter¬ 
rain  data  deck  directly  following  the  start  deck  (10F8.1  format). 
Whether  the  terrain  data  are  on  cards  or  tape,  the  first  terrain 
height  data  point  is  assumed  to  be  located  in  the  southwest  corner 
of  the  grid  and  ordered  such  that  the  elements  in  each  row  are 
read  from  west  to  east,  and  the  rows  are  read  from  south  to  north. 
The  terrain  data  statement  option  exists  for  2-D  problems  in  which 
the  volume  of  terrain  data  is  minimal  and  may  be  "hardwired"  for 
a  specific  terrain  configuration. 


Figure  C.3  Free-Form  Input  Deck  Setup. 
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GLOSSARY  OF  INPUT  VARIABLES 


VARIABLE 

LOCATION 

TYPE 

DEFAULT 

DEFINITION 

ALFA1 

1 

Real 

0.70710678 

Defines  horizontal  transmissivity 
Th  -  1.0/(2.0*ALFA1**2) 

ALFA2 

2 

Real 

0.70710678 

Defines  vertical  transmissivity 

Tv  -  1.0/(2.0*ALFA2**2) 

ATCP 

343 

Real 

Wind  direction  at  top  of  mesh 
(degrees) 

BETA 

3 

Real 

0.1429 

Exponent  of  power  law  profile 

DBUG 

5 

Integer 

0 

Debug  flag 

0:  no  debug  edit 

1:  turn  on  debug  edits 

DMP 

39 

Real 

0.0 

Damping  factor  for  conformal 
simulations  of  highly  stable 
flow;  suggest; 

DMP  -v.05  -*•  .1  for  rv  <  .01 

DX 

6 

Real 

2000.0 

Constant  zone  size  in  x-direction 
(east)  (m) 

DY 

7 

Real 

2000.0 

Constant  zone  size  in  y-direction 
(north)  (m) 

DZ 

8 

Real 

Constant  zone  size  in  z-direction 
(vertical)  (m).  Set  for  real 
space  runs  only 

DZMIN 

9 

Real 

30.0 

Minimum  vertical  zone  size  (m). 
Specified  when  ISIG  *  3 

EPSLN 

10 

Real 

0.1 

min  A<J>.  ..  ■*  conformal  space 

1JK 

min  Div^  real  space 

I2D 

23 

Integer 

0 

Space  dimension  flag 

0:  3-D  space  (x,y,z  or  o) 

1:  2-D  space  (x,y  or  o) 


I 


VARIABLE 

LOCATION 

TYPE 

DEFAULT 

DEFINITION 

rDAY 

11 

Integer 

1 

Day  of  month  being  simulated 

DEBIT 

12 

Integer 

0 

Edit  flag 

0:  no  edit  of  flow  variables 

1:  edit  flew  variables 

2:  generate  line  printer 
contour  plots 

3:  generate  sunnary  edit 
plus  line  printer 
contour  plots 

IHOUR 

13 

Integer 

1 

Hour  of  day  being  simulated 

IMAX 

14 

Integer 

30 

Maxinun  number  of  cells  in  x- 
direction  (east) 

IMN 

15 

Integer 

1 

Minimun  column  in  I -direct ion  to 
be  edited 

IMX 

16 

Integer 

30 

Maximum  column  in  I -direction  to 
be  edited 

IPATH 

17 

Integer 

2 

Geometry  flag 

1 :  real  space 

2 :  conformal  space 

ISIG 

18 

Integer 

3 

Vertical  zoning  flag  (conformal 
space  only) 

1:  log  linear  distribution 

2 :  constant  zoning 

3 :  geometric  progression 

ITERAN 

19 

Integer 

3 

Terrain  data  source  flag 

1 :  read  terrain  data  from  tape 
(unit  8)  [8F10.8  format] 

2:  initialize  terrain  data  from 
data  statement 

3:  read  terrain  data  from  cards 
[10F8. 1  format] 

nvAx  ^ 

40 

Integer 

100 

Maximum  allowable  number  of 
iterations 
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VARIABLE 

LOCATION 

TYPE 

EEFAULT 

DEFINITION 

imm 

20 

Integer 

0 

Special  edit  flag 

0:  no  special  edit 

1:  interpolate  adjusted  field 
to  define  wind  speed  and 
direction  at  specified 
locations  within  the  mesh 

IWIND 

21 

Integer 

1 

Wind  data  type  flag 

1:  surface  wind  data  only, 
Input  wind  speed  and 
direction 

2:  surface  wind  data  only, 
input  wind  components 

3:  surface  wind  data  and  one 
upper  air  station,  input 
wind  speed  and  direction 

4:  surface  wind  data  and  one 
upper  air  station,  input 
wind  components 

IYEAR 

22 

Integer 

1 

Simulation  year 

JMAX 

24 

Integer 

30 

Maximum  nuitoer  of  cells  in 
y-direction  (north) 

JMN 

25 

Integer 

1 

Minimum  column  in  J-direction  to 
be  edited 

JMX 

26 

Integer 

30 

Maximum  column  in  J-direction  to 
be  edited 

KMAX 

27 

Integer 

15 

Maximum  number  of  cells  in  vertical 
direction 

KMN 

28 

Integer 

1 

Mininun  row  to  be  edited 

KMX 

29 

Integer 

15 

Maximun  row  to  be  edited 

yotm 

30 

Integer 

1 

Month  of  simulation 

MPS 


Windspeed  units  flag 

0:  input  is  mi/h 
1:  input  is  m/sec 
2:  input  is  knots 


31 


Integer 


1 


VARIABLE 

LOCATION 

TYPE 

DEFAULT 

DEFINITION 

NDAYS 

32 

Integer 

1 

Number  of  separate  simulations 
to  be  performed 

NPTS 

33 

Integer 

1 

Number  of  wind  observation 
stations 

-1:  initialization  with 
upper  air  data  only 

NSTAT 

291 

Integer 

0 

Number  of  stations  where  wind  speed 
and  wind  direction  are  to  be  com¬ 
puted  for  editing  purposes.  If 
NSTAT  is  negative  4  character 
Hollerith  station  I.D.'s  are  to 
be  entered  via  card  input  (see 
Section  4.1) 

NUA 

345 

Integer 

2 

Number  of  upper  air  reporting 
levels 

QMEG 

34 

Real 

-1.0 

Over  relaxation  parameter  set  to 
-1.0  in  DEFINE,  used  initially 
as  a  flag 

-1.0:  QMEG  is  calculated 
>0  CMBG  was  entered 

through  cards 

ORIOC 

ORIGNY 

35.36 

Real 

0. 0,0.0 

Position  of  grid  origin  (m) 

OfTOP 

342 

Real 

999.0 

Wind  speed  at  top  of  mesh.  If 
QfTCI^999.0  wind  speed  at  top  of 
mesh  is  assured  equal  to  that 
at  top  of  boundary  layer  (m/sec, 
mi/h,  knots) 

ZBL 

344 

Real 

200.0 

Boundary  layer  height  (AGL)  (m) 

ZDIST 

38 

Real 

2*[Raage  of 
Terrain 
Heights] 

Height  of  top  boundary  with  respect 
to  mininun  terrain  altitude  (m) 

UAMSL 

396 

Real  Array 

Array  of  upper  level  heights  (MSL) 
corresponding  to  reported  upper 
level  winds  (m) 

WDR 

41 

Real  Array 

Array  of  observed  surface  wind 
direction  or  V-ccnponents 
(degrees,  m/sec,  mi/h,  or  knots) 

WDRUA 

371 

Real  Array 

Array  of  observed  upper  air  wind 
direction  or  V-ccnponents  read 
from  lowest  level  to  highest  level 
(degrees,  m/sec,  mi/h,  or  knots) 

C-14 


VARIABLE 


LOCATION 


TYPE 


DEFAULT 


DEFINITION 


WSPD 

91 

Real  Array 

Array  of  observed  surface  windspeeds 
or  U-components  (m/sec,  ml/h,  or 
knots) 

WSPDUA 

346 

Real  Array 

Array  of  observed  upper  air  wind 
speeds  or  U-ccmponents  read  from 
lowest  level  to  highest  level 
(m/sec,  mi/h,  or  knots) 

xo 

141 

Real  Array 

Array  of  X-positions  of  observed 
wind  velocities  (m) 

XSTA 

292 

Real  Array 

X- locations  of  each  of  NSTAT 
stations  (m) 

YO 

191 

Real  Array 

Array  of  Y-positions  of  observed 
wind  velocities  (m) 

YSTA 

317 

Real  Array 

Y- locations  of  each  of  NSTAT 
stations  (m) 

ZD 

241 

Real  Array 

Array  of  Z-positions  (AGL)  of 
observed  wind  velocities  (m) 

> 


* 

* 
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SAMPLE  CALCULATION 


To  facilitate  familiarization  with  the  code,  it  is  use¬ 
ful  to  discuss  a  sample  problem.  The  test  case  chosen  is  a 
3-D  Gaussian  Hill  in  the  presence  of  a  westerly  wind  as  shown 
in  Figure  C.4.  Input  conditions  must  include  control  flag 
specifications,  zoning  requirements,  terrain  height  data,  and 
the  initial  velocity  profile. 


rtreaowla*  distance 


Figure  C.4  Simulation  Geometry  for 
Gaussian  Hill  Sample 
Problem. 

The  simulation  was  run  in  the  conformal  mode  with  vari¬ 
able  sigma  spacing  defined  by  a  geometric  progression  assuming 
a  minimum  zone  size  at  the  surface  equal  to  0.013882  and  an 
expansion  factor,  y  =  1.2.  Horizontal  zoning  was  taken  to  be 
750  m  in  both  the  x  and  y  directions.  The  Gaussian  Hill  ter¬ 
rain  height  data  follows  from 
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z(x,y)  -  1000.00  EXP  - 


(C.l) 


u-V> 


2  ♦  (y-yQ)2 

(3000) 2 


where  x  =  y  =  9375.0.  A  total  of  25  zones  was  used  in  both 
o  o 

the  x  and  y  directions,  with  15  zones  specified  in  the  vertical. 
The  velocity  profile  was  defined  by  a  0.14  power-law  up  to  200  m 
and  was  held  constant  and  equal  to  10  m/sec  above  200  m.  The 
top  of  the  mesh  was  taken  to  be  3000  m,  and  a  neutral  atmosphere 
was  assumed. 


Given  these  problem  specifications  and  the  glossary 
of  input  variables,  the  input  deck  structure  may  be  defined. 
Figure  C.5  represents  the  start  deck  and  shows  the  use  of  the 
free-form  variable  format  defined  in  earlier  sections.  The 
values  of  the  terrain  height  array  are  given  by  Eq .  (C.l)  and 
are  tabulated  in  Figure  C.6.  The  overall  deck  structure  is 
shown  in  Figure  C.7. 


•ZO  $241,200.0 
•WSPD  $91,10.0 
•WES  $41,270.0 
•JMN  $25,12 
•JMX  $26,13 
•IMN  $15,1 
•IMX  $16,25 
•KMN  $28 . 1 
•KMX  $29,13 
•  I EDIT  $12,1 
•DX  $6,750.0 
•DY  $7,750.0 
•IMAX  $14,25 
•JMAX  $24,25 
•KMAX  $27,12 
•ZDIST  $38,1500.0 
•DZMIN  $9,50.0 
•OMEG  $34,1. 78 
•EKD  OF  DATA  $ 


Figure  C.5  Input  Deck  Structure. 
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Figure  C.6  Tabulated  Terrain  Data  for  Sample  Problem. 


Figure  C.6  Tabulated  Terrain  Data  for  Sample  Problem  . 
(continued) 
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Using  the  input  deck  structure  shown  in  Figure  C.5  and 
the  terrain  data  from  Figure  C.6,  VARMET  first  calculates  the 
finite  difference  mesh  arrays  (Figure  C.8)  and  the  initial  flow 
profile  (Figure  C.9).  These  data  are  printed  in  subroutines 
SETUP  and  OBSWND,  respectively,  and  are  useful  in  verifying  the 
initial  conditions.  Given  the  initial  wind  field,  the  velocity 
adjustment  phase  is  entered.  Here  the  velocity  potential  field 
is  first  calculated.  Figure  C.10  indicates  the  iteration  by 
iteration  reduction  in  the  maximum  potential  residual.  Once  the 
velocity-potential  field  is  known,  adjustments  to  the  initial 
velocity  field  are  made.  The  final  velocities,  potentials,  and 
divergences  are  then  printed  in  subroutine  EDIT.  An  example  of 
the  line-printer  output  is  shown  in  Figure  C.ll  and  represents 
only  a  few  columns  about  the  hill  peak.  Here  U,V,W  are  the 
velocity  potential  and  DIV  the  local  divergence.  This  edit  is 
arranged  so  that  field  variables  are  edited  from  the  top  of  the 
mesh  down.  The  value  999.0  or  999  indicates  that  an  edit  of  a 
variable  would  be,  not  applicable. 
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MESH  PARAMETERS 


LOCATION  OP  ORIC1I*  (METERS) 

ORIClfX*  «.  OR I CRY*  0.  ORICRZ*  I . 3230E-03 

DX=  7 . 3000E+02  DY*  7.30OOE+O2  ZD 1ST*  3. 0000 E+ 03 

MAXIMUM  HEIGHT  OF  TERRAIN  (METERS)*  1.0000E+O3 


X- ARRAY  (METERS) 


1 

3 • 7300E+02 

2 

1 . 1230E+03 

3 

1 . 8750E+03 

4 

2 . 6250E+03 

3 

33750E+03 

6 

4. 1230E+03 

7 

4. 8750E+03 

8 

3 . 6230E+O3 

9 

6 . 3750E+03 

10 

7. 1250E+03 

11 

7.8750E+03 

12 

8.6250E+03 

13 

9 . 3750E+03 

14 

1.0123E+04 

15 

1 . 0B73E+04 

16 

1 .  1625E+04 

17 

1 . 2375E+04 

18 

1 .3123E+04 

19 

1 . 3873E+04 

20 

1.4625E+04 

21 

1 . 5375 E+04 

22 

1 . 6 123E+04 

23 

1 . 4675 E+04 

24 

1 . 7625E+04 

23 

1 .8373 E+04 

Y- ARRAY  (METERS) 

1 

3 . 7300E+02 

2 

1 .  1230E+03 

3 

1 . 87S0C+03 

4 

2.6250E+03 

3 

3 . 3730E+03 

4 

4. 1250E+03 

7 

4 . 8750E+03 

a 

3 . 6250E+03 

9 

6 . 3750E+03 

10 

7.  1230E+03 

1 1 

7 . 8750E+03 

12 

8. 6250E+03 

13 

9 . 3750E+03 

14 

1 .0125E+04 

15 

1 . 0875E+04 

16 

1 . 1625 E+04 

17 

1 . 2375E+04 

18 

1 .3125E+04 

19 

1 . 3875E+04 

20 

1.4625 E+04 

21 

1 . 5373E+04 

22 

1.6125E+04 

23 

1 . 6875E+04 

24 

1 . 7625E+04 

23 

1  . 8373E+04 

Z- ARRAY  (METERS) 

1 

2. 7326E+03 

2 

2 . 2425E+03 

3 

1 . 8340E+03 

4 

1 . 4937E+03 

5 

1 .2100E+03 

6 

9 . 7363E+02 

7 

7.7666E+02 

8 

6. 125 lE+02 

9 

4 . 7572E+02 

10 

3.  6 173E+02 

1 1 

2 . 6673E+02 

12 

1 . 8737E+02 

13 

1 . 216 1E+02 

14 

6 . 6634E+0 1 

13 

2. 0823E+01 

SICMAF 

ARRAY 

1 

8.6325E-06 

2 

I.7024E-01 

3 

3 . 2677E-0 1 

4 

4.5054E-0I 

S 

5.5369E-01 

6 

6.39G4E-01 

7 

7. 1 127E-01 

8 

7. 7096E-0 1 

9 

8. 2070E-0 1 

10 

0. 6213E-0 1 

11 

8.9670E-01 

12 

9 . 2548E-0 1 

13 

9 . 4047E-O 1 

14 

9 . G946E-0 1 

15 

9.8612E-01 

16 

1 . 0000 E+ 00 

OSIGMAF  ARRAY 

l 

1 .7823E-01 

2 

1 .4853E-01 

3 

1 . 2377E-0 1 

4 

1.0314E-0I 

5 

8.5954E-02 

6 

7. 1623E-02 

7 

5 . 9690E-02 

8 

4.9742E-0C 

9 

4.  1431E-02 

10 

3.43OE-02 

1 1 

2.8786E-02 

12 

2.3988E-02 

13 

1 . 9990E-02 

14 

1 . 6658E-02 

15 

1 . 3882E-02 

16 

I.3882E-02 

SICMAC 

ARRAY 

l 

8.9125E-02 

2 

2.3251E-01 

3 

3 . 8866E-0 1 

4 

5.0212E-01 

5 

5 . 9666E-01 

6 

6 . 7546E-0 1 

7 

7.411  IE-01 

8 

7.9583E-01 

9 

8. 4143E-0 1 

10 

8. 7942E-01 

1  1 

9.  1 109E-01 

12 

9 . 3748E-0 1 

13 

9.39-16E-01 

14 

9.7779E-01 

(5 

9.9306E-01 

DSIGMAC  ARRAY 

1 

1 . 6338E-0 1 

2 

1.633OE-01 

3 

1.3615E-01 

4 

I. 1346 E-Ol 

3 

9.4549E-02 

6 

7. 879  IE-02 

7 

6 . 5659E-02 

a 

5.4716E-02 

9 

4.3597E-02 

10 

3 . 7997E-02 

1 1 

3.  1664E-02 

12 

2.6387E-02 

13 

2. I909E-02 

14 

1 .0324E-02 

15 

1 . 5270E-02 

Figure  C.8  Mesh  Parameters  for  Gaussian  Hill 
Sample  Calculation. 


C-22 


SIGMA  U  V 

8.9I25E-02  9 . 99 13E+00  7. 0296E-07 

2.323  IE-01  9 . 9733E+00  1 . 8900E-06 

3 . 0666E-0 1  9 . 9623E+00  2.8763E-06 

5.0212E-0I  9 . 93 19E+00  3.6963E-06 

5 . 9666E-0 1  9 . 9423E+00  4.3733E-06 

6 . 7346E-0 1  9 . 9353E+00  4.9456E-06 

7.411  IE-01  9.9290E+00  5.4177E-06 

7 . 9323E-0 1  9 . 9233E+00  5.3107E-06 

8. 4 I43E-0 1  9-9194E+00  6.1378E-06 

0.7942E-01  9 .9 I58E+00  6.4102E-06 

9 . 1 I09E-0 1  9 . 9 127E+00  6.6371E-06 

9 . 3740E-0 1  9.9 106E+00  6.7964E-06 

9 . 3946E-0 1  9 . 3272E+00  6.3963E-06 

9 . 7779E-0 1  0.3733E+00  5.8796E-06 

9 . 9306E-9 1  7.283iE*00  4. 996! E-OS 


Figure  C.9  Sample  Problem  Initial 
Flow  Profile. 


ITERATION 

NUMBER* 

1 

MAXIMUM 

RESIDUAL* 

2.9277E+03 

ITERATION 

NUMBER* 

MARI MUM 

RES I DUAL* 

1 .  1301E+03 

ITERATION 

NUMBER* 

3 

MAN I MUM  RESIDUAL* 

7.2749E+02 

ITERATION 

NUMBER* 

4 

MAXI HUM  RESIDUAL* 

5 . 9778E-^02 

ITERATION 

NUMBER* 

3 

MAXIMUM 

RESIDUAL* 

4 . 53 10E+02 

ITERATION 

NUMBER* 

6 

MAXIMUM 

RESIDUAL* 

3.4307E+02 

ITERATION 

NUMBER* 

7 

MAXIMUM 

residual* 

2.7330E-02 

ITERATION 

NUMBER* 

O 

U 

MAXIMUM 

RESIDUAL* 

2.2632E+02 

ITERAT'D’! 

NUMBER* 

9 

MAXIMUM  RESIDUAL* 

I . 96C4E+02 

ITERATION 

NUMBER* 

10 

MAXIMUM 

RESIDUAL* 

l . 6942E^02 

ITERATION 

NUMBER* 

1 1 

MAXIMUM 

RESIDUAL* 

1 . 4340E+02 

ITERATION 

NUMBER* 

12 

MAXIMUM 

RESIDUAL* 

1 . 2 1 E+02 

ITERATION 

NUMBER* 

13 

nAXirnm 

RESIDUAL* 

1 .0I30E-*-02 

ITERATION 

NUMBER* 

14 

MAXIMUM 

RESIDUAL* 

0.  3566E-*-0 1 

ITERATION 

NUMBER* 

13 

fLAxmmi 

RESIDUAL* 

6 . e233E^0 I 

ITERATION 

NUMBER* 

16 

MAXIMUM 

RESIDUAL* 

5 . 532QE+0 1 

ITERATION 

NUMBER* 

• 

• 

17 

• 

• 

MAXIMUM 

RESIDUAL* 

• 

9 

4 . 450GE+0 I 

• 

9 

ITERATION 

• 

NUMBER* 

• 

84 

MAXIMUM 

* 

RESIDUAL* 

• 

2.6I44E-02 

ITERATION 

NUMBER* 

83 

maximum 

RES IDUAL* 

2.3360E-02 

ITERATION 

NUMBER* 

C6 

MAXIMUM 

RESIDUAL* 

2. 1 232E-02 

ITERATION 

NUMBER* 

87 

MAXIMUM 

RESIDUAL* 

1 .9 134E-02 

ITERATION 

NUMBER* 

83 

riAXirnri 

RESIDUAL* 

1 . 7243E-02 

ITERATION 

NUMBER* 

CO 

MAXI HUM 

RESIDUAL* 

1 . 5339E-02 

ITERATION 

NUMBER* 

90 

n.Axirrm 

RESIDUAL* 

1 . 4003E-02 

ITERATION 

NUMBER* 

9  1 

MAXIMUM 

RESIDUAL* 

1 . 26 I9E-02 

ITERATION 

NUMBER* 

92 

MAXIMUM 

RESIDUAL* 

1 .  1372X-02 

ITERATION 

NUMBER* 

on 

MAX  TRIM 

RESIDUAL* 

1 .0242E-02 

ITERATION 

NUMBER* 

94 

MAXIMUM 

RESIDUAL* 

9 . 2352E-03 

ITERATION  CONVERGED 

ITER*  94  RESIDUAL*  9.2352E-03 

Figure  C.10  Maximum  Residual  in 
Potential  Field. 
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Figure  C.ll  Edit  of  Flow  Field  for  Columns  Adjacent  to  Hill  Peak. 


